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ABSTRACT 


A  thermodynamic  system  model  of  a  direct 
injection-type  diesel  engine  has  been  developed. 

The  objective  of  the  modeling  was  to  enable  predic¬ 
tions  of  the  instantaneous  temperature  distribution 
along  and  into  the  combustion  chamber  of  a  diesel 
engine . 

The  present  model  has  a  main  emphasis  on  the 
heat  transfer  computation  with  inclusion  of  new  sub¬ 
models  of  conduction,  radiation  heat  transfer  of  the 
engine,  plus  engine  energy  conversion  processes. 

The  development  of  submodels  of  space-resolved  heat 
release  -and  species  .distribution  in  the  global  model 
was  not  possible  during  the  period  of  the  present 
contract  (01  August,  1984  -  31  July,  1985) . 
Consequently,  some  assumptions  were  employed  for  re¬ 
placing  the  said  submodels  in  the  selected  sample 
computations  based  on  some  of  the  reported  results 
on  heat  release  and  heat  transfer  rates  made  avail¬ 
able  in  literature. 

Some  results  obtained  from  the  present  computa¬ 
tional  study  are  included  in  the  report  and  the  use¬ 
fulness  of  the  new  computer  program  for  the  actual 
engine  design  is  suggested  herein.  In  addition, 
recommendation  for  improvement  of  the  present  compu¬ 
tational  model  is  included. 
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The  recent  adiabatic  engine  activities  have  greatly 
been  enhanced  by  both  excitement  with  the  engine's  superior  per¬ 
formance  and  and  difficulty  of  problems  encountered  in  leading 
the  new  adiabatic  engine  concept  to  an  engineering  product. 
Technical  problems  with  such  stimulation  have  always  been 
favorite  targets  of  many  researchers  in  the  past.  With  juicy  ex¬ 
citement  in  one  hand  and  the  difficulty  in  other  hand,  however, 
the  problems  may  not  be  readily  solved  due  to  their  special 
nature.  For  instance,  the  tribology  problem  of  the  new  engine 
system,  one  of  the  most  detrimental  issues  faced  in  its 
development,  may  be  effectively  solved  only  through  interdis¬ 
ciplinary  collaborations.  A  successful  selection  of  material  and 
its  processing  for  the  construction  of  the  cylinder  liner  needs  a 
careful  analysis  in  various  aspects,  i.e.,  the  surface 
temperature,  the  chemical  composition  of  the  working  environment, 
etc.  In  particular,  when  the  designer's  choice  is  to  employ  a 
composite  material  for  .the  liner,  either  its  experimental  or 
theoretical  analysis  (e.g.,  in  the  above  aspects)  becomes  no 
simple  task. 

One  of  the  more  effective  tools  used  by  the  engine 
designers  under  similar  circumstances  has  been  the  computer 
simulation  of  the  combustion  system.  Such  an  approach  did  not 
provide  the  exact  answer  to  the  problems  but  it  did  greatly  save 
development  costs  in  many  cases.  The  pay-off  of  using  computa¬ 
tional  methods  will  become  greater  when  more  realistic  models  are 
constructed  in  particular  for  the  analysis  of  highly  involved 
problems . 


Most  of  the  previous  engine  analysis  models  developed 
for  a  diesel  engine  are  unable  to  predict  space-resolved  engine 
information,  such  as  the  temperature  distribution  over  the  piston 
crown.  This  was  mainly  due  to  the  difficulties  in  developing 
meaning  ful  models  of  radiation  and  convection  heat  transfer 
processes.  In  addition,  modeling  of  transient  thermal  conduction 
in  engine  components  has  been  severely  limited  in  the  past  com¬ 
pounding  the  said  difficulties.  The  present  work  was  under  taken 
to  facilitate  the  solutions  of  the  issues  delineated  herein. 
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1.  BACKGROUND 


The  high  payoff  development  of  diesel  engine  technology 
has  become  even  more  promising  recently  with  the  introduction  of 
the  innovative  adiabatic  diesel  engine  concept.  U.S.  activities 
in  the  adiabatic  engine  area,  largely  pioneered  by  TACOM,  have 
brought  about  a  significant  technological  transformation  in  the 
diesel  engine. 

A  successful  adiabatic  (the  term,  "low-cooled”  or  "un¬ 
cooled"  may  be  more  suitable,  instead)  engine  is  expected  to  of¬ 
fer  tremendous  advantages  over  the  conventional  diesel  engine. 

The  most  important  operational  characteristics  of  the  new  engine 
may  be  its  uncompromisingly  high  fuel  efficiency,  as  demonstrated 
by  a  low-cooled  diesel  engine  installed  in  a  five-ton  military 
truck  which  achieved  a  brake  specific  fuel  consumption  (bsfc)  of 
0.285  lb/hp-hr  (Bryzik  and  Kamo,  1983).  The  bsfc  is  expected  to 
attain  up  to  0.25  lb/hp-hr  when  the  minimum  friction  ringless 
turbocompound  adiabatic  engine  becomes  available  around  the  year 
1995.  This  new  concept  of  engine  system  certainly  serves,  to  the 
authors*  opinion,  as  a  vitality  injected  into  the  currently 
rather  sluggish  progress  of  efficiency  improvement  in  conven¬ 
tional  diesel  engines:  at  the  beginning  of  their  increasing 
popularity  (in  the  late  1970 ' s) ,  diesel  engine  passenger  cars 
provided  about  29%  higher  fuel  economy  than  similar  cars  with 
gasoline  engines;  the  advantage  of  the  cflesel  passenger  cars,  by 
1983,  decreased  to  only  20%  relative  to  the  gasoline  engine- 
powered  cars  (Wade  and  Jones  1984) .  One  of  the  more  promising 
approaches  to  restoring  the  name  of  "diesel-efficient  engine”  is 
recommended  to  take  the  low-cooled  turbocompound  engine  concept 
in  the  diesel  engine.  The  need  for  improvement  is  equally  recog¬ 
nized  for  large  scale  engines. 

Among  other  comparably  advantageous  characteristics  oi 
the  new  system  are  the  elimination  of  the  cooling  system,  reduced 
size  and  weight,  improved  reliability,  greater  vehicle/engine 
design  freedom,  reduced  noise,  and  greater  fuel  tolerance.  Since 
these  are  important  engine  characteristics  in  meeting  the  propul¬ 
sion  system  goal  of  the  U.S.  Army  vehicle  system,  the  central 
role  of  the  Army  in  the  U.S.  low-cooled  engine  programs  seems  to 
be  entirely  justified.  The  principal  merits  of  the  new  engine 
are  extremely  valuable  in  both  commercial  and  military  applica¬ 
tions.  For  example,  the  use  of  the  uncooled  engine  in  a  combat 
vehicle  will  reduce  the  engine  size  by  50%,  the  cool \ng  system  by 
66%,  and  the  fuel  system  by  40%,  projecting  a  volume  reduction  of 
the  total  propulsion  system  by  39%  (Glance  and  Hunt,  1982) .  In 
addition,  since  up  to  50%  of  the  commercial  and  military  engine 
field  failures  are  related  to  the  cooling  system,  the  absence  of 
a  cooling  system  in  the  engine  will  greatly  reduce  the  vul¬ 
nerability  of  the  vehicle  system.  Furthermore,  since  the  new  en¬ 
gine  offers  a  high  specific  power  density  (hp/ft3) ,  the  system  is 
considered  to  be  a  strong  candidate  for  replacing  the  high  fuel 
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consumption  gas  turbine  presently  used  in  helicopters  (Wilstead, 
1982) .  It  is  expected  the  replacement  will  achieve  more  than  a 
42%  bsfc  reduction  in  existing  helicopters. 

The  low  cooled  diesel  engine  which  is  an  intriguing 
concept  has  created  a  world-wide  challenge  for  engineers,  i.e., 
solving  the  problems  encountered  in  leading  this  concept  to  a 
finished  engineering  product.  As  implied  by  its  name,  the  low 
cooled  engine  has  the  goal  of  controlling  heat  transfer  in  the 
system.  This  may  be  achieved  by  taking  one  of  or  a  combination 
of  the  following:  (1)  proper  development/selection  of  engine 
material,  (2)  creative  design/combination  of  engine  components, 
and  (3)  develop  high-performance  lubrication  for  the  new  engine 
operations.  Since  this  goal  is  presently  being  achieved  by  the 
modification  and/or  redesign  of  the  conventional  direct  injection 
type  engine,  it  is  an  important  requirement  in  the  process  of  its 
development  to  properly  evaluate  the  impact  created  by  the  new 
engine  modification  (here,  the  insulation  and  elimination  of  the 
cooling  system)  on  the  conventional  engine's  function  and 
components.  The  results  of  the  evaluation  are  expected  to  serve 
as  decisive  input  in  the  selection  of  the  component  of  material, 
design,  and  lubrication  for  achieving  a  better  low  cooled  diesel 
engine.  One  of  the  more  useful  means  of  such  evaluation  was  con¬ 
sidered  to  employ  a  computational  method:  an  extensive  engine 
simulation  model  including  the  processes  of  cycle  analysis, 
combustion,  heat  transfer,  etc. 

In  the  construction  of  a  system  analysis  model  of 
diesel  combustion  for  predicting  the  temperature  distribution  in 
the  engine  parts  comprising  the  combustion  chamber,  the  following 
computational  models  that  constitute  the  global  model  are  needed. 

(1)  Engine  cycle  simulation  model, 

(2)  Engine  combustion  model  which  provides  the  dis¬ 
tribution  of  in-cylinder  species, 

(3)  Conduction  heat  transfer  model, 

(4)  Convection  heat  transfer  model, 

(5)  Radiation  heat  transfer  model,  etc. 

At  the  time  when  the  present  work  was  proposed,  there 
were  some  activities  being  carried  out  by  others  in  an  attempt  to 
achieve  a  similar  goal  to  the  present  objectives.  In  spite  of 
such  others'  work,  it  was  decided  to  develop  a  new  system 
analysis  model  in  our  laboratory  because  of  the  following  reason: 
it  was  considered  that  our  new  radiation  heat  transfer  model 
which  had  been  constructed  was  the  most  advanced  of  the  kind;  and 
we  have  been  carrying  out  the  proposed  study  by  developing  other 
submodels  in  the  above  list.  Further  the  proposed  study  was 
closely  related  with  and  comprises  our  ongoing  engine  combustion 
research  at  Rutgers,  the  objective  of  which  was  to  obtain  the 
relationship  of  radiation  heat  transfer  to  in-cylinder  species 
distribution  in  a  diesel  engine. 
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2.  OBJECTIVE 


The  objective  of  the  present  study  was  to  construct  a 
computational  system  analysis  model  of  a  direct  injection-type 
diesel  engine.  The  model  was  to  enable  prediction  of  the  tem¬ 
perature  distribution  into  and  along  the  combustion  chamber  of  a 
diesel  engine  made  of  composite  material.  The  original  goal  of 
the  model  use  was  to  implement  the  analysis  by  employing  an  IBM 
PC  so  that  some  handy  analysis  could  be  carried  out  was  needed. 


3.  METHODOLOGY 


Among  the  submodels  of  engine  processes  mentioned  in 
the  previous  page  are  two  submodels  that  were  being  developed  by 
the  time  of  its  proposal  in  our  laboratory.  They  were  submodels 
of  the  engine  cycle  simulation  and  the  radiation  heat  transfer 
process.  In  particular,  the  radiation  model  had  actively  been 
constructed  on  an  ARO-supported  sister  project.  They  were  in  a 
usable  condition  but  needed  some  work  of  debugging  and  adjustment 
for  being  included  in  the  present  global  model  of  engine  system 
analysis  to  be  performed  during  the  project  period. 

As  an  overall  evaluation  of  the  methodology  taken  in 
the  present  study,  although  we  hoped  to  construct  a  global  model 
of  engine  system  analysis  as  perfect  as  possible  during  the 
project  period,  it  turns  out  that  we  have  taken  an  approach  or 
achieving  bilateral  objectives  in  its  development.  One  was  a 
short-term  goal  to  construct  a  new  thermal  conduction  model, 
which  has  become  the  most  interesting  and  intriguing  result  among 
those  obtained  in  the  present  project.  The  key  issue  of  the  con¬ 
duction  model  is  a  new  numerical  scheme,  called  "the  control 
volume  finite  element  method"  developed  on  the  present 
sponsorship.  Because  of  our  emphasis  on  its  new  numerical  scheme 
due  to  its  potentially  extensive  and  great  usefulness  in  various 
ways,  the  global  model  construction  was  made  with  inclusion  of 
more  assumptions  than  we  originally  planned  for.  Consequently, 
the  construction  of  an  more  in-depth  model  of  the  combustion 
process  was  placed  as  a  long-term  goal  to  be  implemented  as  a 
continuing  activity  in  our  laboratory.  In  describing  the  results 
obtained  on  the  present  sponsorship,  the  report  is  prepared  in 
two  parts: 

(1)  A  detailed  discussion  of  the  new  conduction  model 
along  with  a  brief  description  of  other  submodels,  and 

(2)  Analysis  of  computational  results  obtained  by  using 
the  present  global  model  in  predicting  the  temperature  distribu¬ 
tion  in  a  piston  of  a  direct  injection-type  diesel  engine. 

Note  that  an  extensive  manuscript,  which  will  contain 
more  details,  is  being  prepared  for  the  purpose  of  reporting  in  a 
technical  society. 


4.  DESCRIPTION  OF  THE  SYSTEM  ANALYSIS  MODEL 

In  order  to  construct  a  system  analysis  model  of  a 
diesel  engine  which  enables  the  predictions  of  temperature  dis¬ 
tribution  in  the  engine  components  of  the  combustion  chamber  and 
of  the  engine  performance,  several  submodels  need  to  be 
developed. 

Combustion  Process.  The  engine  combustion  process  is 
in  need  of  proper  modeling.  Its  modeling  may  be  achieved  by 
modeling  of  various  processes  including  the  gas  exchange  (intake 
and  exhaust)  process,  the  fuel  injection/atomization,  mixing  and 
vaporization,  the  gas  motions,  the  chemical  kinetics,  etc.  The 
construction  of  such  a  model  is  important  and,  on  the  other  hand, 
is  difficult  to  achieve:  it  is  important  to  construct  a  model 
because,  in  a  computation  of  heat  transfer  in  the  cylinder,  for 
example,  such  a  model  will  provide  the  source  of  heat,  i.e.,  the 
temperature  and  species  distributions;  it  is  difficult  to  model 
the  process  due  to  so  many  unknowns  and  complexity  of  the  various 
processes,  i.e.,  modeling  "unsteady  heterogeneous  reacting  tur¬ 
bulence  flows."  Even  with  some  success  in  modeling  of  gas 
motions,  its  coupling  with  the  energy  equation  has  been  a  dif¬ 
ficult  task  to  overcome.  In  spite  of  the  difficulties,  some 
success  in  its  modeling  is  reported  in  recent  literature  (Amsden, 
et  al.,  1985).  Again,  the  inclusion  of  such  a  model  in  our 
global  model  would  have  been  highly  desirable.  Note,  that, 
however,  the  construction  of  such  a  model  or  an  endeavor  of  com¬ 
bining  an  existing  combustion  model  to  a  heat  transfer  model  is  a 
lengthy  task  requiring  in-depth  experience  in  modeling  of  the 
kind.  Under  the  given  limitations  in  time  and  resource  faced  in 
the  proposed  project,  an  approximation  model  of  combustion 
process  was  to  be  developed  for  the  global  model.  This  choice 
included  in  the  original  plan  was  abandoned  in  order  to  con¬ 
centrate  our  effort  for  a  thorough  study  of  heat  transfer  rather 
than  diverting  a  part  of  the  resource  to  construct  a  model  which 
is  sure  to  be  obsolete  even  upon  completion  of  its  construction. 
Instead,  assumptions  were  employed,  which  will  be  removed  or  min¬ 
imized  in  a  next  study,  for  the  global  model  as  explained  later. 

Once  the  heat  source,  which  may  be  characterized  by  the 
instantaneous  distributions  of  temperature  and  species  is  known, 
a  fairly  extensive  heat  transfer  computation  can  be  made  by  using 
the  present  global  model.  The  convection  model  may  become 
straight  forward,  at  least  in  theory,  upon  the  completion  of  the 
combustion  modeling:  its  validity  will  be  as  good  as  that  of  the 
combustion  model.  The  radiation  heat  transfer  modeling  is  one  of 
the  main  parts  for  the  global  model.  In  addition,  the  conduction 
heat  transfer  model  may  be  placed  as  the  main  program  of  the  heat 
transfer  model  to  predict  the  temperature  distribution  over  the 
combustion  chamber.  Computational  results  from  the  other  sub¬ 
models  may  be  introduced  as  the  boundary  condition  of  the  main 
program . 


RADIATION  HEAT  TRANSFER  MODEL.  Since  the  fraction  of 
heat  transfer  through  the  radiation  process  out  of  the  total  heat 
loss  in  a  typical  direct  injection  type  diesel  engine  may  be  30% 
or  more  (Ebersole,  et  al.,  1963),  it  is  necessary  to  include  the 
radiation  heat  transfer  analysis  in  a  system  analysis  of  the 
engine.  In  particular,  since  the  radiation  heat  transfer  process 
is  highly  direction  and  spectrum  dependent  in  nature,  its  proper 
analysis  is  needed  in  predicting  the  temperature  distribution 
over  a  diesel  engine  combustion  chamber.  In  view  of  its 
importance,  a  computational  and  experimental  research  project  is 
being  carried  out  at  Rutgers  under  a  separate  sponsorship  of  ARO. 
A  brief  discussion  is  made  on  the  modeling  of  radiation  heat 
transfer  in  diesel  combustion. 

In  implementation  of  a  better  analysis  of  radiation 
heat  transfer  in  diesel  combustion,  several  difficult  problems 
are  faced:  Most  serious  problems  to  be  solved  in  modeling  the 
process  may  be  listed  as  follows: 

(1)  A  detailed  description  of  distributions  of  combus¬ 
tion  products  and  temperature  in  terms  of  a  suitable  coordinate 
system  centered  at  individual  locations  over  the  combustion  cham¬ 
ber  wall.  This  information  is  necessary  for  finding  the  optical 
properties  in  the  hemispherical  volume  faced  by  each  location, 
which  are  used  in  solving  the  equation  of  radiation  through  in¬ 
dividual  optical  paths  with  respect  to  the  location.  • 

(2)  A  viable  computation  of  spectral  volume  absorptance 
of  mixture  of  combustion  products.  This  is  required  for  process¬ 
ing  the  distributions  mentioned  in  (2)  to  find  the  optical 
properties  to  be  used  in  the  governing  equation  of  radiation  to 
be  mentioned  later. 

(3)  A  more  accurate  computational  method  for  integrat¬ 
ing  the  equation  of  radiation  heat  transfer  along  individual  op¬ 
tical  paths  in  the  hemispherical  volume  as  stated  in  (1) . 

In  order  to  overcome  the  abovementioned  difficulties 
the  following  methods  were  developed/employed:  Mentioning  the 
problem  stated  in  (1) ,  the  species  distribution  was  computed  by 
using  a  new  coordinate  transformation  technique  developed  in  our 
laboratory  (Chang  and  Rhee,  1985;  a  copy  of  its  manuscript  is  in¬ 
cluded  in  Appendix-I) .  The  key  of  the  new  method  is  to  find, 
upon  the  species  distribution  (temperature  distribution,  as  well) 
given  in  terms  of  a  coordinate  with  respect  to  the  injection 
nozzle,  the  species  distributions  along  individual  optical  paths 
centered  at  any  locations  over  the  combustion  chamber  wall. 
Speaking  of  the  issue  explained  in  (2) ,  the  volume  absorptance  of 
combustion  products  were  calculated  for  the  range  of  radiation  by 
using  the  Rayleigh  limitation  expression  and  a  semiempirical  band 
model  in  finding  soot  radiation  and  gas  radiation,  respectively. 
This  is  discussed  in  more  detail  in  the  abovementioned 


manuscript.  A  better  computational  technique  for  solving  the 
governing  equation  of  radiation  heat  transfer  was  sought  by  using 
the  Gaussian  integration  method  (Chang  and  Rhee,  being  prepared 
for  publication).  More  discussion  on  this  is  made  in  Appendix-I. 
Note  that  the  results  from  the  radiation  heat  transfer  computa¬ 
tion  are  needed  for  finding  boundary  conditions  to  be  introduced 
in  the  conduction  heat  transfer  computation. 

CONDUCTION  HEAT  TRANSFER.  Since  this  subject  carries  a 
more  weight  than  others  in  the  present  report,  a  separate  chapter 
is  assigned  for  its  extensive  discussion  as  follows. 


5.  THREE-DIMENSIONAL  CONDUCTION  HEAT  TRANSFER  MODELING  FOR 
CYLINDRICAL  SHAPE  BODIES 

The  problem  of  finding  the  heat  transfer  and  tempera¬ 
ture  distributions  in  the  combustion  chamber  has  long  been 
studied.  It  is  important  to  have  the  distributions  because  the 
engine  manufacturers  face  strong  needs  of  a  proper  thermal  load¬ 
ing  analysis  of  engine  parts  composing  the  combustion  chamber.  A 
number  of  approaches  were  taken  to  measure  and  predict  the  tem¬ 
perature  distribution  in  such  complex  geometry  components.  Dis¬ 
cussing  the  experimental  method,  in  1964,  Whitehouse  et  al. 
measured  the  variation  in  a  piston  crown.  He  found  that  the 
fluctuation  in  temperature  at  the  'surface  was  rapidly  reduced  to 
steady  states  within  the  depth  of  a  few  milimeters  of  the 
surface.  Woschini  and  his  coworkers  (1967;  1979;  1979)  measured 
the  temperature  of  a  piston  at  20  discrete  points  and  evaluate 
the  entire  steady  state  temperature  distribution  by  using  both 
the  relaxation  method  and  electrolytic  tank  analysis.  Regarding 
the  computational  analysis,  the  finite  element  method  has  often 
been  employed  to  obtain  the  temperature  distributions.  Saursten 
(1979)  and  Saugerud  and  Sandsmark  (1979)  assumed  a  steady  steady 
axisymmetric  thermal  boundary  condition  and  calculated  the  dis¬ 
tribution  in  a  piston.  They  also  attempted  a  comparison  of  a 
two-dimensional  computation  and  results  measured  from  a  real 
three  dimensional  case.  Recently,  Morel  et  al.  (1984)  used  a 
one-dimensional  transient  conduction  code  to  analyze  the  tempera¬ 
ture  fluctuation  in  a  thin  layer  of  the  combustion  chamber  and 
used  the  finite  element  method  to  calculate  the  steady  state  tem¬ 
perature  distribution  inside  the  chamber  wall.  Completing  a 
literature  survey,  it  is  found  that  there  has  not  been  any 
suitable  three-dimensional  transient  conduction  analysis  con¬ 
ducted  for  the  engine  components.  For  example,  the  number  of  en¬ 
gine  cycles  from  the  start  of  an  engine  until  various  critical 
places  in  its  combustion  chamber  reach  steady  state  temperature 
has  not  been  properly  evalyuated.  Recognizing  that  such  a  piece 
of  information  is  of  great  value  for  accessing  the  engine'  ther¬ 
mal  loading,  a  new  extensive  was  developed  by  taking  the  follow¬ 
ing  approaches. 


[•] 


(1)  A  new  control  volume  finite  element  method  that 
holds  the  merits  of  both  the  finite  element  and  finite  difference 
methods  was  developed. 

(2)  A  successive  differential  method  in  time  domain  was 
employed  to  obtain  more  accurate  and  stable  solutions. 

(3)  A  new  computer  aided  automatic  mesh  generation 
program  which  expedites  the  implementation  of  the  model  was 
developed . 

The  new  conduction  model  was  developed  to  be  placed  as 
the  main  program  in  the  present  global  analysis  of  heat  transfer 
in  an  engine.  Results  computed  from  other  submodels  are  intro¬ 
duced  into  the  main  program  as  boundary  conditions.  In  the 
present  study,  the  time-  and  space-resolved  temperature  distribu¬ 
tion  in  a  diesel  engine  piston  was  calculated  by  using  the  bound¬ 
ary  conditions,  as  mentioned  in  Chapter  4.  Further  details  of 
the  processes  taken  to  obtain  this  are  explained  in  the 
following. 


A.  DESCRIPTION  OF  THE  NEW  CONDUCTION  ANALYSIS  METHOD 

The  problem  given  here  is  to  find  the  instantaneous 
temperature  distribution  in  the  engine  parts  of  a  direct  injec¬ 
tion  type  diesel  engine  when  boundary  conditions  are  known.  The 
geometry  of  the  engine  components  is  highly  irregular  in  shape 
and  due  to  the  stratification  nature  of  combustion  reactants  and 
products,  the  boundary  conditions  vary  with  location  and  time. 
Such  boundary  conditions  are  sought  in  other  submodels  as  pre¬ 
viously  mentioned. 

The  general  heat  conduction  problem  with  no  heat  gener¬ 
ation  in  a  region,  R ,  enclosed  by  a  boundary  is  governed  by  the 
following  differential  equation. 

dT  _ 

pc  —  +7  •  q  -  0  (1) 

dt 

where,  q  =  -k^T,  p  is  the  mass  density  if  material,  c 
is  the  specific  heat  and  k  is  the  thermal  conductivity.  The  de¬ 
pendent  variable,  T  is  continuous  in  R ,  and  quantities  p,  c  and 
k  are  uniform  in  each  subregion  of  R .  There  may  be  interfaces 
at  which  p ,  c  and  k  are  discontinuous  as  long  as  the  heat  flux 
and  temperature  are  continuous  over  the  interface.  What  has  been 
don  here  is  to  develop  a  new  numerical  technique  written  in  a 
computer  program  to  analyze  an  unsteady  three-dimensional  conduc¬ 
tion  process  as  expressed  by  Eq.  (1) .  in  engine  components  with 
cylindrical  geometry.  In  order  to  solve  the  stated  problem  the 
following  tasks  were  carried  out. 
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(a)  Automatic  Mesh  Generation, 

(b)  Control  Volume  Approach, 

(c)  Matrix  Formation  and  Manipulation,  and 

(d)  Solution  of  the  Discritization  Equation. 


Automatic  Mesh  Generation  in  Two-dimensional  Plane. 

One  of  the  first  steps  to  be  taken  in  numerical  analysis  of  heat 
conduction  may  be  the  generation  of  meshes  in  the  computation 
domain  with  appropriate  nodal  points  representing  the  individual 
nodal  elements.  Often,  this  process  is  tedious  and  prone  to  gen¬ 
erate  errors:  development  of  a  computer  aided  mesh  generation 
technique  was  highly  desirable  to  use  in  our  work  so  as  to  handle 
elements  with  complex  geometry.  The  technique  sought  was 
required  to  meet  the  following  criteria: 

(1)  Adequate  boundary  description  indicating  if 
connected  singulary  or  multiply; 

(2)  Means  to  identify  material  properties  in  each  node; 

(3)  Acceptable  resolution  and  accuracy; 

(4)  Easy  bookkeeping  of  coordinate  of  nodes;  and 

(5)  Easy  graphical  presentation  by  using  a  plotter. 

A  new  mesh  generation  technique  has  been  developed  in  a 
direction  to  meet  the  above  requirement  as  much  as  possible.  The 
program  was  written  in  Basic  language  to  be  used  in  an  IBM  PC 
connected  to  an *HP  7475  plotter  to  handle  elements  with  two- 
dimensional  region.  Since  our  new  numerical  scheme  utilizes 
nodal  elements  with  triangular  geometry,  the  mesh  generation  was 
made  to  meet  this  purpose. 

The  essence  of  the  automatic  mesh  generation  scheme  is 
the  use  of  piecewise  connection  of  successive  lines  from  nodal 
point  to  nodal  point  in  order  to  triangulated  the  computation 
domain.  Consider  a  typical  geometry  of  substance  of  heat  trans¬ 
fer  which  is  to  be  divided  into  nodal  elements  over  the  entire 
domain  as  shown  in  Fig  1.  The  domain  geometry  may  be  either  ir¬ 
regular  or  regular.  The  domain  will  be  divided  into  quadrangular 
elements  with  more  likely  uniform  sizes  with  its  points  along  the 
boundary  contour  successively  connected.  Subsequently,  the  quad¬ 
rangle  is  divided  into  two  triangular  elements,  in  either  way  as 
shown  in  elements  P  and  Q  of  Fig.  1. 

One  of  the  main  input  data  is  a  set  of  sufficient  num¬ 
ber  of  coordinates  showing  the  points  to  describe  the  boundary 
contour  of  the  domain.  The  only  requirement  is  the  number  of 
points  to  be  specified  have  to  equal  each  other  in  both  sides  and 
both  ends,  respectively;  the  number  in  one  side  does  not  need  to 
be  the  same  as  that  in  one  end,  however.  In  determining  the  ele¬ 
ments  of  quadrangle  with  relatively  uniform  sizes  in  the  domain 
the  following  method  is  employed.  By  taking  an  arbitrary  node 
(i,j),  one  may  define  the  following  new  numbers. 


rli  =  <xl,j  “  xl,l>/<xl,n+l  -  xl,l> 

r2j  "  <xm+l# “  xm+l,l)/(xm+l,n+l  "  xm+l,l) 

sli  "  ^il  "  “  ym+1,1^ 

s2j  "  (Yi,m+1  "  Ym+i,n+l)/(yi,n+l  _ym+l,n+l) 


where,  m  is  the  total  number  of  points  in  one  side  and  n  is  the 
number  of  points  along  one  end.  The  determination  of  nodal 
points  ,yji)  within  the  contour  of  the  substance  are  deter¬ 
mined  by  usingJthe  following  formulas: 


(m  + 

1-  i)/(m  +  r2j 

(i 

-  l))/m 

Ij 

(n  + 

1  -j)/(n  +  s2j 

(j 

-  l))/n 

xij  ~  rij)  xil  +  rij  xin 
Yij  -  (1  -  si:j)  ylfj  +  si;j  yij. 

During  the  process  of  developing  the  computer-aided 
mesh  generation  by  taking,  the  abgve  method,  several  refinements 
were  made  for  its  uses  in  analysis  of  various  geometries.  Among 
the  significant  refinements  is  a  special  caution  for  the  case 
with  right-angled  corners  of  the  domain  to  be  meshed.  Without 
this,  in  such  a  case,  multiply  defined  points  or  ill-conditioned 
elements  (e.g.,  triangles  with  zero  area)  were  encountered.  An 
example  of  results  from  the  automatic  mesh  generation  is  shown  in 
Fig.  2.  The  number  in  each  nodal  point  indicates  the  corner 
characteristics  which  determines  a  pertinent  numerical  scheme  as 
explained  later.  Without  mentioning  the  details,  note  that  the 
characteristic  numbers  are  determined  in  the  program  which 
governs  the  energy  conservation  equation  to  be  explained  later. 


Modal  Elements  in  Cylindrical  Shape  Elements.  For  the 
determination  of  nodal  elements  in  cylindrical  bodies,  the  r-z 
plane  of  the  component  was  divided  into  triangular  elements  by 
using  the  abovementioned  two-dimensional  method.  Upon  the  deter¬ 
mination  of  meshes  in  r-z  plane,  the  component  was  divided  in  the 
circumferential  direction  as  shown  in  Fig.  3. 


B.  CONTROL  VOLUME  FINITE  ELEMENT  METHOD. 


Since  the  exact  solution  of  Eq.  (1),  in  general,  is  not 
readily  available  if  the  geometry  of  the  domain  becomes  complex, 


numerical  methods  are  employed,  e.g.,  the  finite  difference 
method,  the  finite  element  method  and  the  control  volume  method. 
Although  the  logic  and  algebra  of  the  finite  difference  method 
are  simple,  its  discrete  node  elements  of  which  geometry  is 
limited  to  square  geometry  often  lack  in  resolution,  e.g.,  the 
incomplete  meshes  near  the  sharp  corner  boundaries.  The  dif¬ 
ferential  equation,  Eq.  (1) ,  for  such  incomplete  meshes  have  to 
be  approximated  for  every  individual  case  making  its  computer 
program  of  general  use.  The  finite  element  method,  on  the  other 
hand,  can  deal  with  irregular  geometries,  since  it  employs  trian¬ 
gular  elements  opposed  to  the  square  elements  in  the  finite  dif¬ 
ference  method.  Mathematically,  the  finite  element  method  is 
based  on  variational  formulation  or  Gallichin  weighted  integral 
formulation  (Shin,  1984) .  The  mathematical  formulation,  however, 
is  rather  involved  and  makes  the  computer  program  much  more  com¬ 
plex  than  that  of  the  finite  difference  method.  Further,  its 
physical  interpretation  is  not  as  clear  as  its  counterpart  in  the 
finite  difference  method.  This  shortcoming  often  outweights  its 
advantages  in  particular  when  used  in  two-dimensional  analysis: 
its  advantage  may  be  recognized  in  analysis  of  three-dimensional 
problems. 

In  the  present  work,  a  control  volume  formulation  is 
employed.  It  combines  the  advantages  of  both  the  finite  dif¬ 
ference  and  finite  element  methods  in  dealing  with  the  domains 
with  irregular  geometry  in  a  simple  manner.  The  control  volume 
approach  is  based  on  the  conservation  of  thermal  energy.  The 
method  is  based  on  the  numerical  balance  of  thermal  energy  over 
small  subelements  within  individual  nodal  elements,  consequently 
over  the  entire  computation  domain.  This  is  achieved  by  starting 
from  an  integral  conservation  statement:  again,  the  statement  ap¬ 
plies  locally  and  globally. 

An  integral  equation  can  be  obtained  based  on  Eq.  (1) 
over  a  domain  as 

■Ir  f  PCTdfl  +  /  q-ndo  =0  (2) 

3t  «  o 

where,  <?  is  the  region  being  considered  and  o  is  the  boundary  of 
the  region.  From  the  integral  equation,  we  can  derive  the  dis¬ 
cretization  equation  for  each  element.  The  basic  procedure  of 
deriving  the  discretization  equation  of  the  present  method  are, 

(1)  discretize  the  calculation  domain  into  elements 
of  certain  geometric  shape, 

(2)  use  the  element  based  interpretation  function, 

(3)  derive  the  discretization  equations  by  using 
the  control  volume  method,  and 

(4)  Compilation  of  coefficients  in  the  discritization 
for  all  the  elements. 


Two-dimensional  Discretization  and  its  Equations.  The 
discretization  of  the  calculation  domain  in  a  two-dimensional 
plane  is  made  as  follows.  As  mentioned  previously,  the  domain  is 
divided  into  triangular  elements.  Consider  a  collection  of  the 
elements  in  the  plane  as  shown  in  Fig.  4-(a).  The  element  again 
divided  into  three-node  elements  by  joining  any  of  the  midpoints 
of  its  side  as  shown  in  Fig.  4-(b).  A  typical  triangular  element 
with  its  subelements  will,  then,  become  one  as  shown  in  Fig.  5. 
Upon  the  determination  of  the  triangular  elements  and  their 
subelements,  the  integral  conservation  statement  is  applied  to 
each  of  the  three  subelements.  When  applied,  the  integral  equa¬ 
tion  for  each  subelement  becomes 


-r-  /  pcTdfl  +  /  q.ndo  +  J  q.ndo 

Ot  n 

+( similar  contributions  from  other  subelements  associated 

with  node  l)+( boundary  conditions,  if  applicable )=0  * 

•|r  /  pcTdfl  +  /  q.ndo 

«, 

+( similar  contributions  from  other  subelements  associated 
with  node  2)+(boundary  conditions,  if  applicable)=0 

/  "pcTdfl  +  /  q.ndo 

3 

+( similar  contributions  from  other  subelements  associated 
with  node.  3 )  +  ( boundary  conditions,  if  applicable)=0 

The  integrands  in  Eq.  (3)  are  approximated  based  on 
averaged  thermal  properties,  e.g.,  kj4  «  2  k^  k-s  /  (k^  +  k^). 
The  above  equations,  therefore,  are  rewritten  as  follows.  J 


a  UAi 

at  '  pcTdfl  =  picA5E 


a 


a  Uia 

at  '  pcTdfi  “ 

W3 


/  q.ndo  = 


T  -  T 
i  A  j 


15 


UMNIHiP  Jl  V1  l  W.mf  }  ITJ  .■  1  >'*  '!•  1  V 
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In  these  equations,  A^-s  represents  the  vertical  area  between  the 
subvolume  i  and  j;  D^-sJdenotes  for  the  distance  between  the 
subelements  i  and  j  as  defined  in  the  following. 


DRij  =  C(Xi  -  Xj)2  +  (yi  -  Yj)2]1/2. 


Introducing  these  approximation  terms  into  the  conser- 
*  vation  equation  to  individual  nodal  points,  the  integral  conser- 

?  vation  equation  of  the  corresponding  triangular  element  becomes, 

for  example,  with  respect  to  node  1, 
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(5)  . 


Regarding  the  elements  composing  the  boundary  of  the 
domain,  the  domain  is  declared  to  hold  a  specific  temperature  of 
heat  flux  replacing  a  term  in  the  above  equation.  The  above 
equation  can  be  written,  in  a  straightforward  manner  in  the  com¬ 
puter  program,  for  the  entire  elements  in  the  domain  with  inclu¬ 
sion  of  appropriate  boundary  conditions.  What  has  to  be  done, 
completing  this  step,  is  to  solve  the  simultaneous  unsteady  equa¬ 
tions  for  finding  temperature  values  at  individual  nodal  points, 
which  will  be  explained  later. 


Three-Dimensional  Discretization  and  Discritization 
Equations.  Although  an  extension  of  the  present  two-di»ensional 
analysis  may  be  made  for  solutions  of  three-dimensional  domain  by 
using  the  tetrahedral  element.  Upon  a  brief  evaluation  of  this 
approach,  it  was  found  that  there  are  a  number  of  problems  to  be 
overcome  in  making  the  ideal  of  practical  use,  which  is  presently 
being  studied  in  our  laboratory.  On  the  other  hand,  since  most 
of  the  engine  elements  comprising  the  combustion  chamber  are 
axisymmetric,  a  new  method  has  been  developed  for  analysis  of 
such  objects,  which  is  considerably  simple  to  implement. 


The  basic  idea  of  the  present  new  method  to  analyze 
three-dimensional  objects  is  as  follows:  first  of  all,  the  r-z 


fchiiVii  3m  ifcnlWln  Ii 
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pane  of  the  cylindrical  coordinate  centered  at  the  axis  of  a 
cylindrical  three-dimensional  domain  is  discntized  by  using  the 
previously  mentioned  new  numerical  scheme  developed  for  the  two- 
dimensional  domain,  and  next  the  circumferential  direction  of  the 
domain  is  discretized  by  appropriate  increments.  Two  kinds  of 
typical  building  blocks  of  nodal  elements  for  this  method  are 
shown  in  Fig.  6.  The  elements  were  divided  into  six  subelements 
in  the  corresponding  nodal  points  and  the  discretization  equa¬ 
tions  of  energy  conservation  are  written.  The  elemental  equa¬ 
tions  for  those  subvolumes  are  listed  in  the  following. 


•r—  /  pcTdfl  +  /  q.ndo  +  /  q.ndo  +  /  q.ndo 
ot  n 

*1  ®1,  °1« 

+  ( similar  contributions  from  other  subelements  associated 

with  node  l)+(boundary  conditions,  if  applicable)=0 


/  pcTdfl  +  /  q.ndo  +  /  q.ndo 

+( similar  contributions  from  other  subelements  associated 
with  node  2 )+( boundary  conditions,  if  applicable )=0 

0  4  •¥  -4 

/  pcTdfl  +  /  q.ndo  +  /  q.ndo 

a  n 

3  VJ1  3  * 

+( similar  contributions  from  other  subelements  associated 
with  node  3 )+( boundary  conditions,  if  applicable )=0 


-r-  /  pcTdfl  +  /  q.ndo  +  /  q.ndo  +  /  q.ndo 

<?  o  o  o 

+( similar  contributions  from  other  subelements  associated 
with  node  4)+(boundary  conditions,  if  applicable)=0 

3  .  -*  -» 

/  pcTdfl  +  /  q.ndo  +  /  q.ndo 

A,  otl  oi4 

.  +( similar  contributions  from  other  subelements  associated 
with  node  5)+(boundary  conditions,  if  applicable)=0 


0  4  4  4  4 

•st-  /  pcTdt?  +  /  q.ndo  +  /  q.ndo 

Ot.  - 

M  a  o 

•  *  3  t* 

-►(similar  contributions  from  other  subelements  associated 
with  node  6 )+( boundary  conditions,  if  applicable )=0 


(6). 


Using  a  similar  approximation  as  that  in  the  previously 
mentioned  two-dimensional  analysis,  one  can  write  the  integral 
terms  in  the  above  equation  in  the  following  concise  form: 
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ij  ID 


(i, j=l,2---6) 


(7) 


Note  that  this  equation  is  applied  to  tie  entire  ele¬ 
ments  with  appropriate  boundary  conditions.  In  Eq.  (7)  ,  A^-s  is 
the  area  normal  to  the  direction  through  which  heat  flow  from 
subcontrol  volume  i  to  j  and  DRjj  is  the  distance  between  the 
two.  Consequently,  a  linear  ordinary  differential  equation  can 
be  established  for  individual  subelements  associate  with  the  cor' 
responding  nodal  point.  It  is  further  found  that  the  final  form 
of  the  discretization  equation  for  temperature  at  nodal  point  i, 

,  becomes , 


dT± 

T7—  =  EaT  +  b  +  b  t 
dt  n  n  n  io  ii 


where  the  summation  is  made  over  all  nodes  pertinent  to  subele¬ 
ments  associated  with  node  i. 


Solution  Algorithm.  In  solving  simultaneous  unsteady 
ordinary  differential  equations,  expressed  as  a  form  of  Eq.  (8) , 
the  equations  are  reduced  to  matrix  form  as  follows, 


A  u  +  B[u  +ut] 

D  D 


u(tc )  =  u0 


(9)  • 


Although  there  are  several  methods  to  solve  such 
equations,  one  way  of  increasing  accuracy  and  stability  in  time 
marching  is  to  employ  "high  order  approximation"  by  using  the 
successive  differentiation  method.  Without  showing  the  details 
of  its  derivation,  the  result  of  the  present  successive  differen 
tiation  method  is  shown  in  the  following. 


U  =  A  u+3(u  +U  t) 
D  D 


u  =  A  u  +B  u  =A ( A  u)+A  B ( u  +ut )+B  u 
D  D  D  o 


«M» _ i  l  _  i _ i  _ 3  _  _  2  _  _  _  i  _ i 

u  =Au  =A  U  +A3U  =A  u-rA  B  t  u  +U  t)+ABu 
d  d  b  b 


(10)  . 


;(i 


M  -i - 


_  i  -i  -l  —  t 
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b  b  b 
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u(t)=u0  +  (Au0+Bub)t+  ,|2  — -u0t+  + 


-i-l-  -i-2- 
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6.  COMPUTER  PROGRAM  FOR  THE  CONTROL  VOLUME  FINI'il  ELEMENT  METHOD 
BY  USING  THE  SUCCESSIVE  DIFFERENTIATION 


A  new  computer  program  has  been  developed  utilizing  the 
new  concept  of  the  present  control  volume  finite  element  method 
by  using  a  numerical  scheme  of  the  successive  differentiation. 

The  program  was  written  in  FORTRAN  language  and  a  name  was  given 
as  "HEATRAN."  The  HEATRAN  program  of  the  heat  conduction 
analysis  along  with  its  flow  chart  is  included  in  Appendix-II, 
where  the  entire  global  model  of  heat  transfer  analysis  is 
included. 


7.  PREDICTION  OF  TEMPERATURE  DISTRIBUTION  IN  A  DIESEL 
ENGINE  PISTON 


The  present  global  model  which  includes  several 
submodels,  e.g.,  a  conduction  model  and  radiation  heat  transfer 
model,  was  used  for  predicting  the  temperature  distribution  in  a 
diesel  engine  piston.  The  following  assumptions  were  employed  in 
the  computation: 

(1)  The  instantaneous  space  averaged  heat  flux  through 
the  piston  top  surface  is  given  as  shown  in  Fig.  7. 

(2)  Of  the  heat  flux  history  in  Fig.  7,  the  distribu¬ 
tion  of  instantaneous  heat  flux  during  the  period,  350  -  400 
crank  angle  (CA) ,  however,  is  space-resolved  as  follows:  the  dis¬ 
tribution  in  r-direction  is  given  in  the  form  as  shown  in  Fig.  8, 
and  the  distribution  in  circumferential  direction  is  of  normal 
distribution.  The  integrated  total  heat  flux  at  each  engine 
crank  angle  given  in  the  abovementioned  way  is  equal  to  that  in 
corresponding  engine  crank  angle  shown  in  Fig.  7. 

(3)  The  total  heat  flux  can  be  dissolved  into  portions 
of  those  through  radiation  and  convection  processes  in  either 
ways  as  shown  in  Fig.  8  or  Fig.  9. 
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(4)  The  piston  is  made  of  aluminum. 


(5)  The  cut-view  of  r-z  plane  of  the  piston  is  shown  in 
Fig.  2  and  it  is  axisymmetric. 

(6)  The  side  of  the  piston  is  insulated.  The  bottom 
surface  of  the  piston  is  maintained  at  400  deg  K.  The  piston  was 
initially  at  400  K. 

(7)  The  engine  is  of  direct-injection  type  and  has  four 
identical  spray  plumes  with  distributions  of  radiatively  par¬ 
ticipating  species,  F,  expressed  in  Eg.  (8)  in  Appendix-I.  The 
engine  has  minimum  swirl  motions. 

The  reasons  for  employing  the  assumptions  may  be  ex¬ 
plained  in  the  following.  The  instantaneous  heat  flux  through 
the  piston  top  surface  in  diesel  combustion  is  not  properly 
known,  although  some  measurements  and  predictions  have  been 
reported,  e.g.,  those  by  Whitehouse  (1970).  Some  active  study  is 
carried  out  in  our  laboratory  for  helping  improve  our  understand¬ 
ing  of  the  process.  Regarding  the  fractions  of  heat  fluxes  con¬ 
tributing  either  convection  or  radiation  process,  the  basic  un¬ 
resolved  issue  is  the  unavailability  of  in-cylinder  species 
distribution.  Since  there  are  two  different  kinds  of  radiation 
heat  transfer  flux  distribution  considered  in  literature  (Chang, 
Yang  and  Rhee,  1985;  Chang,  Kobayashi  and  Kamimoto,  1985)  the 
present  computation  was  carried  out  for  both  cases  (Fig.  8  and 
9) .  The  rest  of  the  assumptions  can  clearly  be  replaced  by  more 
dependable  input  boundary  conditions  for  the  main  program.  More 
realistic  computations  are  being  planned  in  our  laboratory. 

For  reference,  note  that  the  computing  time  required 
for  a  semi-steady  state  computation  starting  from  the  initial 
condition  was  about  five  hours  when  an  AS/9000  computer  located 
at  Rutgers  was  used;  the  results  shown  in  this  report  are  those 
obtained  after  completing  computation  of  about  700  cycles  of 
the  engine  operation  at  a  speed  of  1,000  rpm. 


Instantaneous  Temperature  at  Various  Locations  on  the 
Top  of  Piston.  The  temperature  history  computed  based  on  the 
heat  flux  input  shown  in  Fig.  8  is  shown  for  various  locations  at 
the  piston  top  (Figs.  10,  11  and  12).  Note  that  the  computed 
temperature_results  are  shown  in  a  nondimens ional  value,  T, 
defined  as  T  =  (T  -  400)  /  400.  The  small  section  of  figure  in 
Fig.  10  and  some  of  the  other  figures  herein  presented  indicates 
the  locations  in  the  piston  top  of  the  corresponding  temperature 
history.  While  the  peak  input  heat  flux  is  given  at  r/R  =  0.5  - 
0.8,  the  highest  temperature  level  was  found  to  occur  at  similar 
locations,  which  may  be  expected.  When  histories  of  both  heat 
flux  and  temperature  are  compared,  it  is  found  that  there  is  a 
delay  in  rise  of  temperature  by  about  50  degrees  of  CA  from  the 
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rise  in  the  input  heat  flux.  Such  a  lag  in  temperature  rise  is  a 
typical  example  of  the  transient  conduction  problems.  Further 
analysis  is  needed  on  this  problem  because  of  its  importance  in 
evaluating  the  thermal  loading  on  the  piston. 


Surface  Temperature  Distribution  over  the  Piston  Top 
Surface.  Figure  11  shows  the  space-resolved  temperature  dis¬ 
tribution  plotted  for  engine  crank  angle  of  400.  This  is  a  time 
before  the  peak  temperature  is  reached  on  the  piston  top  so  that 
the  surface  distribution  continues  to  change  as  seen  in  Fig.  10. 
After  brief  look  at  Figs.  10  an  11,  one  can  readily  find  that 
there  is  relatively  small  variation  of  temperature  level  in  cir¬ 
cumferential  direction  although  the  heat  flux  input  was  given  in 
a  normal  distribution  in  the  direction  (see  assumption- (2)  shown 
in  the  previous  page) .  This  may  be  due  to  both  the  high  thermal 
conductivity  of  piston  piston  material  (aluminum)  and  the  rapid 
heat  flows  contributed  by  two  spray  plumes  into  the  space  with 
low  heat  input  located  between  them  (consider  the  valley  between 
two  normal  distributions  of  heat  input) .  The  temperature  dis¬ 
tribution  along  the  direction  of  axis  of  the  plume,  i.e.,  the 
radial  direction  of  the  piston  is  shown  in  another  form  at  engine 
crank  angle  of  400  as  shown  in  Fig.  12.  The  results  in  this 
figure  show  that  the  peak  in  temperature  does  not  correspond  to 
the  peak  found  (see  point  A  in  Fig.  10)  in  the  later  engine  crank 
angle.  Note  however  when  a  computation  was  made  for  the  heat  in¬ 
put  given  as  Fig.  9,  which  has  higher  radiation  heat  flux  near 
the  injection  nozzle,  the  peak  of  the  temperature  was  found  cor¬ 
respondingly  near  the  nozzle  as  well  (see  Fig.  13) . 


Temperature  Distribution  in  an  r-z  Pane  of  the  Piston. 

Without  any  particular  reason,  again  the  results  com¬ 
puted  for  engine  crank  angle  of  400  were  plotted  in  order  to  show 
the  temperature  distribution  in  r-z  plane  of  the  piston.  The 
temperature  distribution  near  the  rings  and  piston  bowl  curves 
are  of  high  detail.  Computational  results  were  plotted  for  the 
heat  input  given  in  Figs.  8  and  9,  for  comparison  purpose,  as 
shown  in  Figs.  14  and  15,  respectively.  Note  that  such  informa¬ 
tion  can  be  readily  evaluated  for  the  pistons  having  the  ceramic 
piston  cap  presently  used  in  the  uncooled  diesel  engine.  It  is 
further  noted  that  when  more  realistic  heat  flux  input  becomes 
available  for  the  locations  near  the  piston  rings,  some  useful 
design  information  may  be  derived  by  using  the  present  model. 
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SUMMARY  AND  RECOMMENDATIONS 


A  preliminary  thermodynamic  system  analysis  model  of  a 
direct  injection- type  diesel  engine  has  been  constructed.  Among 
the  significant  achievements  made  in  the  present  project  are  (1) 
development  of  a  new  numerical  method  of  transient  conduction 
heat  transfer  analysis,  (2)  an  introduction  of  the  successive 
differentiation  method  in  time-marching  in  search  of  its  numeri¬ 
cal  solution,  (3)  construction  of  a  computer  aided  mesh  gener¬ 
ation  program,  (4)  an  analysis  of  temperature  distribution  into 
and  over  a  diesel  engine  piston,  etc.  Further  improvement  of  the 
present  model  is  being  carried  out  and  at  the  same  time  more  in- 
depth  analysis  of  the  temperature  distribution  in  other  engine 
component  of  the  engine  combustion  chamber  is  carried  out. 

Since  the  computation  of  the  temperature  distribution 
in  engine  components  of  the  chamber  is  made  based  on  the  results 
from  the  radiation  heat  transfer  and  species  distribution  in  the 
combustion  chamber,  the  accuracy  of  temperature  analysis  is  ob¬ 
viously  dependent  upon  their  results.  The  logical  consequence, 
therefore,  is  that  there  is  an  absolute  need  in  improving  the 
computation  to  obtain  those  results. 

It  is  recognized  that  due  to  the  complexity,  and  scale 
of  the  undertaken  project  (consider  that  the  computer  time  for  a 
simple  analysis  by  using  our  program  in  an  AS/9000  was  about  five 
hours)  a  long-term  support  for  it  is  needed.  It  is  felt  that 
completion  of  the  present  contract  simply  serves  as  a  starting 
point  for  both  further  analysis  of  the  temperature  distribution 
and  construction  of  a  better  global  heat  transfer  model  of  a 
diesel  engine. 

Since  the  results  shown  in  the  report  may  be  a  mere 
demonstration  to  show  what  the  present  model  can  do  for  achieving 
a  better  analysis  of  an  engine's  heat  transfer,  any  concluding 
technical  remark  on  its  process  may  be  somewhat  premature. 
Nevertheless,  we  believe  that  our  new  model  of  engine  heat  trans¬ 
fer  analysis  can  profitably  be  used  even  in  the  present  form  al¬ 
though  there  is  a  lot  of  room  for  improvement. 


ACKNOWLEDGMENT 


The  present  work  was  proposed  to  conduct  on  the  strong 
encouragement  by  Dr.  Walt  Bryzik  of  U.S.  Army  Tank  Automotive 
Command  (TACOM) .  The  entire  course  of  the  study  was  closely 
monitored  by  Dr.  David  M.  Mann,  U.S.  Army  Research  Office.  Mr. 
Ernest  Schwarz  of  TACOM  served  as  scientific  liason  repre¬ 
sentative  on  the  project.  Their  interests  and  suggestions  on  the 
work  are  greatly  appreciated. 


10.  LIST  OF  PUBLICATIONS  AND  TECHNICAL  REPORTS 

Presently,  two  manuscripts  are  being  prepared  for  pub¬ 
lishing  in  professional  journals. 

11.  PARTICIPATING  SCIENTIFIC  PERSONEL 

The  following  individuals  have  contributed  in  carrying 
out  the  present  project: 

Professor  K.  T.  Rhee,  Principal  Investigator 

Dr.  S.  L.  Chang,  Researcher,  and 

Mr.  X.  L.  Yang,  Graduate  Student. 

Mr.  Yang  is  scheduled  to  have  his  final  examination  for  his  M.S. 
degree  which  has  been  pursed  on  a  partial  finalcial  support  from 
the  present  contract. 


12.  REFERENCES 


Amsden,  A.  A.,  Ramshaw,  J.  D. ,  Cloutman,  L.  D.  and 
O'Rourke,  P.  J.,  "KIVA:  A  Computer  Program  for  Two-  and  Three- 
Dimensional  Fluid  Flows  with  Chemical  Reactions  and  Fuel  Spray," 
LA-10245-MS  and  LA-10534-MS,  Los  Alamos  National  Laboratory,  Sep 
tember  1985. 

Bryzik,  W.  and  Kamo,  R. ,  "TACOM/Cummins  Adiabatic  En¬ 
gine  Program,"  SAE  Paper-830314,  1983 

Chang,  S.  L.  and  Rhee,  K.  T. ,  "Coordinate  Transforma¬ 
tion  Method  for  Radiation  Heat  Transfer  Prediction  in  Soot  Laden 
Combustion  Products,"  being  reviewed  for  publication,  1985. 

Chang,  S.  L. ,  Yang,  X.  L.  and  Rhee,  K.  T. ,  "A 
Parametric  Analysis  of  Radiation  Heat  Transfer  in  Direct  Injec¬ 
tion  Diesel  Combustion,"  Int'l  Symposium  on  Diagnostics  and 
Modeling  of  Combustion  in  Reciprocating  Engines,  Tokyo,  Japan, 
September,  1985. 

Chang,  Y.  J. ,  Kobayashi,  H. ,  Matsuzawa,  K,  and 
Kamimoto,  T.,  "A  Photographic  Study  of  Soot  Formation  and  Combus 
tion  in  a  Diesel  Flame  with  a  Rapid  Compression  Machine,"  Int'l 
Symposium  on  Diagnostics  and  Modeling  of  Combustion  in 
Reciprocating  Engines,  Tokyo,  Japan,  September  1985. 


"The 


Ebersole,  G.  D. ;  Myers,  P.  S.  and  Uyehara,  0.  A., 
Radiant  and  Convective  Components  of  Diesel  Engine  Heat 
Transfer,"  SAE  Paper-701C,  1963. 

Glance,  P.  C.  and  Munt,  R. ,  "Engine  Research 
Development,  and  Acquisition  Strategy  for  U.S.  Army  Ground 
Vehicles,"  Engine/Fuel  Workshop  sponsored  by  U.S.  ARO  and  U.S. 
AMEDRC,  San  Antonio,  Texas,  December  1982. 

Sarsten,  A. ,  "Survey  of  Theoretical  and  Experimental 
Evaluation  of  Thermo  Loading  of  Diesel  Engines  in  Norway,"  SAE 
Paper-790819,  1979. 

Saugerud,  0.  T.  and  Sandsmark,  N. ,  "Strength  Analysis 
of  Thermally  Loaded  Engine  Components  by  Twi-dimensional  and 
Three-dimensional  Finite  Element  Models,"  SAE  Paper-790820, 
1979. 


Shin,  T.  M. ,  Numerical  Heat  Transfer.  Hemisphere  Pub¬ 
lishing  Corp.,  1984 

Wade,  W.  R.  and  Jones,  C.  M. ,  "Current  and  Future  Light 
Duty  Diesel  Engine  and  Their  Fuels,"  SAE  Paper-84015,  1984 

Whitehouse,  N.  D. ,  "Heat  Transfer  in  a  Quiescent  Cham¬ 
ber  Diesel  Engine,"  Proc.  Instn.  Mechanical  Engineers,  185,  963, 
1970. 

Wilsted,  H.  D.,  "Preliminary  Survey  of  Possible  Use  of 
the  Compound  Adiabatic  Diesel  Engine  for  Helicopters,"  SAE  Paper- 
820432,  1982 

Woschini,  G.,  "A  Universally  Applicable  Equation  for 
Instantaneous  Heat  Transfer  in  the  Internal  Combustion  Engine," 
SAE  Paper-670931,  1967. 

Woschini,  G.,  "Prediction  of  Thermal  Loading  of  Super¬ 
charged  Diesel  Engine,"  SAE  Paper-790821,  1979. 

Woschini,  G. ,  Fieger,  J. ,  "Determination  of  Local  Heat 
Transfer  Coefficients  at  the  Piston  of  a  High  Speed  Diesel  Engine 
for  Evaluation  of  Measured  Temperature  Distribution,"  SAE  Paper- 
790834,  1979. 


Locations  along  r  -  Direction 


Point  A 


Instantaneous  Temperature  at  Various  Locations  over  the  Piston  Top 


at  Engine  Crank  Angle  of  400. 


Temperature  Variation  along  r  -  Direction  on  the  Piston  Top 
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Fig.  13  Temperature  Variation  along  r  -  Direction  on  the  Piston  Top 
at  Engine  Crank  Angle  of  400. 
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COORDINATE  TRANSFORMATION  METHOD  FOR 
RADIATION  HEAT  TRANSFER  PREDICTION  IN 
SOOT  LADEN  COMBUSTION  PRODUCTS 

by 

S.  L.  Chang  and  K.  T.  Rhee 
Department  of  Mechanical  and  Aerospace  Engineering 
Rutgers,  The  State  University  of  New  Jersey 
New  Brunswick,  New  Jersey  08903 

INTRODUCTION 

Computational  modeling  of  radiative  heat  transfer  has  recently  become  a 
more  important  means  for  the  analysis  of  thermal  loading  of  combustion  chambers 
with  luminous  flames.  This  seems  to  arise  due  to  both  the  renewed  recognition 
of  need  for  an  extensive  radiation  heat  transfer  analysis  to  achieve  a  better 
design  of  many  practical  systems  (e.g.,  furnace,  turbine  combustor  and  diesel 
engines)  and  promising  new  contributions  from  recent  studies  enabling  more 

m 

compreiiensi ve  modeling.  In  the  modeling,  however,  several  difficult  problems 
are  faced,  mainly  due  to  the  directionality  and  spectral  nature  of  radiative 
processes.  Most  serious  problems  to  be  solved  in  improving  its  modeling  may  be 
listed  as  follows: 

(1)  A  viable  computation  of  spectral  volume  absorptance  of  a  mixture  of 
combustion  products; 

(2)  A  more  accurate  computational  method  for  integrating  the  equation  of 
radiation  heat  transfer  along  individual  optical  paths  in  the  hemi spherical 
volume  faced  by  each  location  of  the  reactor  wall  (hereafter  called  a  detector) 
and  over  the  wave  length  range  of  thermal  radiation;  and 

(3)  A  detailed  description  of  in-reactor  distributions  of  combustion 
products  and  temperature  in  terms  of  a  suitable  coordinate  system  centered  at 
an  individual  detector  and  its  coupling  with  the  above  radiative  transfer- 


equation. 


In  the  present  paper,  being  mainly  concerned  with  the  last  item  of  the 
list,  a  brief  discussion  is  given  for  others  in  order  to  review  the  overall 


problem. 


SPECTRAL  VOLUME  ABSORPTANCE  OF  SOOT/GAS  MIXTURE 

When  the  spectral  volume  absorptance  of  combustion  products  are  calculated 
for  the  range  of  thermal  radiation,  the  respective  absorptances  due  to  soot 
cloud  and  gas  mixture  are  computed  as  well  as  their  summation.  The  most  widely 
accepted  method  of  computing  the  absorptance  of  soot  r  ,  by  using  either 
Rayleigh  limit  expression  [1]  for  median  size  of  particles  smaller  than 
0.1  micrometer,  or  Mie  theory  [2]  for  particles  with  greater  sizes  where  it  is 
required  to  know  the  size  distribution  to  take  into  account  scattering.  Since 
practical  combustion  systems,  in  general,  produce  small-size  soots,  the  Rayleigh 
expression  is  only  considered  here.  The  expression  is  independent  of  the  size 
of  the  particles,  inversely  proportional  to  the  wavelength  \,  and  proportional 
to  the  soot  volume  fraction,  i.e., 


36  n  k(n/\)  f 


2  2  2  4  2 

[n  (1-k  )  +  Z]c  +  4  nV 


(1) 


The  optical  constants,  n,  k,  the  real  and  imaginary  parts  of  the  refractive 
indices,  were  measured  from  the  intensity  of  polarizing  light  reflected  from  the 
soot  particles  and  the  measurements  were  fitted  to  the  dispersion  equation  [3,4] 
It  was  reported  that  *•  may  slightly  increase  with  wavelength  and 
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insignificantly  varies  with  kind  of  soot  and  temperature.  For  computation  of 
emissions  from  gaseous  mixtures,  the  empirical  wide  band  model  proposed  by 
Edwards  and  Bal arkri shnan  [5]  has  often  been  used  [6-8].  The  significance  of 
the  method  is  to  find,  instead  of  finding  gray  body  absorptance  of  gas  volumes. 


3. 


the  width  of  individual  emission  bands,  A^,  compared  to  blackbody  spectrum  to 
approximate  the  strength  of  original  emission  spectrum.  This  unique  technique 
facilitates  an  evaluation  of  the  spectral  dependence  of  gas  radiation.  When 
emission  bands  overlap  each  other  at  a  single  spectrum,  a  separate  correction  is 
given  [5].  With  regard  to  the  summation  of  both  emissions  due  to  soots  and 
gases,  even  though  any  method  has  not  fully  been  verified,  several  techniques 
have  been  used  in  the  past:  for  an  assumed  homogeneous  and  gray  body  emission, 
by  adding  the  gas  emittance  times  the  soot  transmittance  to  soot  emittance  [9]; 
by  adding  the  soot  and  gas  emittances  e$  and  e^,  to  find  total  emittance,  e^, 
i.e.,  Ej.  =  1— (1-c)  ( 1- e$  )  [10];  by  introducing  a  new  approximation  using  a 
pentagamma  functional  form  [11],  etc.  Unlike  those  for  computation  of  gray  body 
emittance  of  gas/soot  mixture  volume,  the  summation  of  spectral  emissivity  of 
such  mixtures  needs  special  consideration  in  order  to  evaluate  the  spectral 
dependence  of  radiation  heat  transfer.  A  new  method  of  summing  spectral  volume 
absorptances  has  been  used  by  the  authors  to  find  the  spectral  emittance  of 
combustion  products  comprising  both  soot  and  gases  for  their  recent  studies  [12]. 
In  the  computation  of  the  spectral  volume  absorptance  of  a  soot/gas  mixture, 
the  spectral  volume  absorptance  of  gaseous  mixture,  ^  is  expressed,  for 
convenience  of  computation,  in  a  form 

r\.g  =  CA  •  (2) 

Here,  C  was  found  from  the  relationship,  A  =  /  1  -Ex p ( — C JL/V )d\  upon  the  computa- 

A  A  A 

tion  of  A  for  each  emission  band  by  using  the  Edwards'  wide  band  non-gray  model. 

A 

On  the  other  hand,  the  volume  absorptance  due  to  soot  suspending  in  the  mixture, 
r  ,  (Eq.  1)  is  used  to  find  C'  from  the  expression  r  =  C'/\.  The  spectral 
absorptance  of  soot/gas  mixture,  r  ,  is  then  computed  from  the  following: 

A 

rx  =  (C+C')A  (3) 


This  method  for  quantifying  spectral  dependence  of  radiative  processes  of 
a  combustion  gas  mixture  with  soot  suspension  is  basically  similar  to  one  of 
those  listed  above  for  gray  body  radiation  analysis.  That  is,  a  summation  for 
individual  spectral  emissivity  of  such  mixtures  by  using  the  new  method,  without 
exceeding  blackbody  emissivity  at  a  spectrum  upon  the  addition  of  both  emissions 
is  quite  equivalent  to  the  following  relationship: 

£t.X  =  ^“'g.xW-'s.x)  <4> 

The  validity  of  the  new  summation  method,  however  similar  to  one  often  employed 
for  gray  body  analysis,  is  not  fully  accepted  at  present,  so  that  a  future 
systematic  experimental  evaluation  of  summation  is  needed.  Further,  since  the 
calculation  of  emittance  as  a  thermal  equilibrium  mixture  is  made  by  computing 
the  spectral  volume  absorptance  instead  of  emittance  as  explained  above,  the 
assumption  of  Kirchhoff's  law  to  -hold  for  it  in  a  system  is  more  justified  than 
a  similar  assumption  made  for  gray  body  computations. 


EQUATION  OF  RADIATION  HEAT  TRANSFER  EQUATION 

The  equation  of  radiation  heat  transfer  in  radiatively  parti ci pati ng  fluids 
confined  in  a  volume  is  well  known.  The  governing  equation  of  monochromati c 
radiation  heat  transfer  is  sought  for  the  processes  in  combustors  containing 
plumes  of  small  size  soot  laden  flames  as  shown  in  Fig.  1  where  pertinent 
details  are  shown,  e.g.,  geometric  information  and  spectral  reflectivity  of  wall 
surfaces,  B^,  etc.  Referring  to  Fig.  1,  the  governing  equation  of  local 
spectral  radiation  intensity,  I  (r),  is  written  as 


f  0  r\(r' 

I  M  =  /  - - 


e  (r1)  Exp(-/  tr  (r")dr")dr‘  , 


where,  e. 


,  r  is  the  spectral  volume  absorptance  at  (r,e,c), 


k  [ Exp( hc/XKT)-l ] 

c  =  2.998  x  1010  cm/sec,  h  =  6.625  x  10  ^  erg-sec,  and  k  =  1.380  x  10  ^  erg/K. 

The  heat  flux  along  an  optical  path  of  a  solid  angle  (e,^)  in  the  thermal  radiation 


range  may  be,  then,  computed  from 


■/„  Ix(0> 


cos  e  dx 


The  main  difficulty  in  finding  the  solution  of  Eq.  (5)  resides  in  the  fact  that 
k  is  a  function  of  both  species  concentrations  and  temperature  and  that  both 
variables  are  a  function  of  the  location  along  the  individual  optical  path,  r.  The 
present  paper  considers  a  new  technique  to  describe  the  species  distribution 
along  individual  optical  paths  and  its  coupling  with  Eq.  (5).  When  such 
distribution  details  become  available,  a  more  realistic  solution  of  Eq.  (6)  is 
accessible,  as  reported  elsewhere  [13]. 

The  surface  properties  are  combined  with  the  radiation  heat  fluxes  incident 
on  individual  locations  of  the  combustor  wall.  The  extinction  coefficient, 

■c  =  /it  •  dL,  along  each  optical  path  found  from  the  computation  as  explained 
above  and  the  spectral  reflectivities,  B,  , ,  are  included  in  the  expression 
for  the  radiosity  of  surface  1,  B. , 

K  i  9  dS 


(I,  +  B,n  cos  en  1.(0))  +  ba1  cos  e.,  e'  (l0  +  b,  .  cos  1,(0)) 


i  -  Ba  bx2  cos  9l  cos  e2  e  x 

Upon  implementation  of  computations  by  coupling  Eq.  (5)  with  Eq.  (7),  the  net  heat 
transferred  through  a  chosen  spectrum  at  a  detector  can  be  calculated. 


IN-COMBUSTOR  SPECIES  DISTRIBUTION 


j  Since  an  analysis  of  space-  and  spectrum-resol ved  radiation  heat  transfer 

in  a  combustor,  as  considered  here,  requires  detailed  distributions  of  in-reactor 
species  and  temperature,  search  of  a  useful  technique  for  this  is  clearly 
|  warranted.  Unfortunately,  at  present  neither  convenient  means  for  obtaining 

in-combustor  data  sufficient  for  use  in  radiation  heat  transfer  analysis  has 
been  available.  The  distributions  may  be  obtained  by  either  experimental  or 
j  theoretical  method.  Some  activity  is  underway  in  the  direction  of  obtaining 

such  information  in  the  authors'  laboratory.  When  the  distribution  becomes 
available,  it  is  most  convenient  to  express  It  in  a  coordinate  system  centered 
I  at  the  fuel  injection  nozzle. 

Although,  as  explained  above,  the  species  distribution  in  plumes  of 
combustors  equipped  with  the  fuel  injection  nozzle  is  not  readily  available,  a 
|  new  approach  explained  in  the  following  may  shed  light  in  the  search  for  a 

viable  solution  to  the  problem.  It  is  found  from  some  experimental  and 
theoretical  results  in  literature  that  the  species  distribution  in 
non-ax i symmetric  plume  of  such  systems  [14,15]  may  conveniently  be  described  by 
the  following  equation  in  a  cylindrical  coordinate: 

F  =  Fq  exp(-ap-b0^-cz^)  (8) 

where  F  may  represent  either  fuel/air  ratio,  or  soot  concentration,  or  CC^  or 
H^O  concentrations;  F  ,  a,  b  and  c  are  constants  to  be  separately  determined  by 
either  experimental  or  theoretical  means  for  individual  combustion  units.  This 
expression  of  plumes,  however,  deemed  rather  simple  and  probably  less  realistic, 
is  considered  to  serve  a  suitable  method  for  some  parametric  analysis  of  the 
confronting  problem.  Among  the  advantages  of  this  approach  is  that  the  deformation 

i 

of  plume,  e.g.  in  spiral  direction,  due  to  swirl  motion  while  the  fuel  spray  is 

I 


formed,  may  be  reasonably  described  for  both  transient  and  steady  combustion 
systems,  i.e,  by  including  it  in  <j>-di  recti  on,  as  explained  later. 

The  next  task  to  be  implemented  for  the  final  solution  may  be  a  proper 
coupling  of,  with  Eq.(5),  the  details  of  in-combustor  species  distribution  that 
may  be  obtained  from  the  above  method  or  others.  Several  techniques  have  been 
employed  in  finding  a  solution  of  the  equation  of  radiation  heat  transfer  { E q . 
(5),  namely,  the  Monte-Carlo  method  [10,16],  the  zonal  method  [17,18],  the  use 
of  geometric  factors  [6,19],  etc.  Those  methods  did  not  rigorously  implement 
the  solution  of  Eq.  (5),  but  sought  for  approximated  solution  and  gray-body 
computation.  Among  the  main  reasons  for  this  is  that  it  is  difficult  to  find 
species  distribution  along  individual  directions  of  integration  in  the 
hemi spheri cal  volume  faced  by  the  location  where  heat  transfer  is  considered. 
The  difficulty  is  compounded  when  the  distribution  varies  with  time  as  in 
transient-flow  combustors. 

COORDINATE  TRANSFORMATION 

In  view  of  the  discussion  above  on  the  overall  problem  in  finding  space- 
and  spectral -resol ved  radiation  heat  transfer  in  combustors,  it  is  highly 
desirable  to  find  a  more  versatile  method  enabling  us  to  find  a  better  solution 
of  Eq.  (5).  A  new  technique  for  the  goal  is  presented  in  the  following.  The 
basic  idea  of  the  technique  is  to  employ  coordinate  transformations  for  finding 
the  species  distribution  along  the  individual  optical  direction,  r,  (Fig.  1)  in 
an  equational  form  by  using  the  species  distribution  given  in  a  coordinate  with 
respect  to  the  fuel  injector.  When  such  an  equation  in  a  suitable  coordinate 
centered  at  detectors  on  the  combustor  wall  becomes  available,  i.e.  along 
individual  optical  direction  with  respect  to  the  detector,  a  better  computation 
of  the  space-resolved  radiation  heat  transfer  can  be  attempted. 


As  explained  earlier,  it  was  found  to  be  convenient  to  express,  in  a 
cylindrical  coordinate,  the  spray  plumes  of  which  geometry  and  others  vary  with 
time.  The  following  discussion  is  concerned  with  description  of  the  same  plumes 
in  a  new  coordinate  system  with  respect  to  a  detector  of  interest,  consequently 
to  make  available  the  species  distribution  along  individual  optical  paths 
centered  at  the  detector.  This  goal  is  achieved  by  using  a  new  coordinate 
transformation  technique  involving  two  main  sets  of  coordinate  systems,  i.e.,  a 
cylindrical  coordinate  centered  at  the  injection  nozzle  hole,  0,  for  describing 
the  spray  plumes  (p,4>,z)  and  a  spherical  coordinate  (r,9,c)  for  expressing  the 
same  plumes  with  respect  to  the  detector  located  at  0.  Two  transformation 
techniques  are  presented:  (1)  point-to-point  transformation;  and  (2)  distribution 
equation-to-distribution  equation  transformati on .  The  former  is  considered  to 
be  particularly  useful  when  some  discrete  point  data  are  available  in  a  coordinate 
with  respect  to  the  injection  nozzle  for  expressing  them  in  a  new  coordinate 
centered  at  a  detector.  The  latter  is  developed  for  transforming  a  functional 
expression  of  the  spray  Dlume  in  a  nozzle  coordinate,  e.g.  Eq.  (8),  to  find  an 
equational  spray  description  in  a  new  detector  coordinate. 

( 1 )  Point-to-point  Transformations 

MUTUALLY  PARALLEL  COORDINATES  -  When  in-combustor  species  data  is 
experimentally  acquired  [20],  It  is  always  likely  that  they  become  available 
in  a  set  of  discrete  data  points  as  opposed  to  in  a  form  of  continuous  variation 
in  species  concentration  of  plumes.  When  one  uses  such  discrete  data  points  in 
computing  radiation  characteri sties  in  a  direction  of  optical  path  with  respect 
to  a  detector,  it  is  necessary  to  know  a  species  distribution  along  individual 
optical  paths  for  solving  Eq.  (5).  For  such  needs,  a  new  approach  of  using 
point-to-point  coordinate  transformat  ion  is  introduced.  To  perform  a  coordinate 
transformation  between  these  systems,  two  more  intermediate  coordinate  systems. 


X  and  Y,  are  introduced  for  the  convenience  of  derivation.  i he  axes  of  the 
coordinates  under  consideration  are  either  parallel  or  perpendicul ar  to  each 
other  in  order  to  avoid  the  complexity  of  the  derivation  although  cases  with 
the  mutually  non-parallel  coordinates  are  discussed  later.  The  interim 
coordinate  systems  are  listed  in  the  following  and  shown  in  Fig.  2-(A). 


a. 

Spray  plume  coordinate  (cylindrical). 

C 

(p» 

♦  .z) 

b. 

X  coordinate, 

X 

(xl 

» x  2  >  x  2 ) 

c. 

Y  coordinate. 

Y 

(yl 

’^2’^  3^ 

d. 

Detection  coordinate  (spherical) 

S 

(r. 

e.O 

Identifying  the  location  of  a  detector  at  x^d»  x2d>  x3d,  a  series  of 
transformations  are  carried  out  to  obtain 

2  2  2  2 
p  =  r  sin  e  +  2  r  pd  sine  cos  (<!>d+0  +  Pd  > 

<t>  =  tan_1[(pd  s i n<t>d  -  r  sine  sin^)/(pd  cos  (td  +  r  sine  cose)], 

Z  =  Zd  -  r  cose  ,  (9) 

2  2  2  -1  -1 
where,  pd  =  Xj  +  x2  >  <t>d  =  tan  (x^x-^).  arid  z d  =  x3-  Note  that  4>  =  tan  (x^/xj) 

if  x1  ^  0  ,  <j>  =  tan”1(x2/x1)  +  *  if  X}  <  0  and  x2>  0  and  <r>  =  tan  1(x2A1)  -  w  if 

X|  <  0  and  x2  <  0. 

MUTUALLY  UNPARALLEL  COORDINATES  -  However,  since  in  many  practical  systems, 
e.g.,  diesel  engines,  the  fuel  is  injected  into  the  combustion  chamber  at  an  angle 
e^,  as  shown  in  Figs.  1  and  2-(B),  this  has  to  be  taken  into  consideration  in 
the  coordination  transformation.  In  order  to  meet  this  need,  an  additional 
coordinate  system  having  an  angle  of  e.n  with  respect  to  the  cylinder 
coordinate  (X|»x2>x2)  is  introduced.  The  new  coordinate  (u^,u9,u3) 
is  not  parallel  to  the  detector  coordinate  (y^y^y^).  Their  mutual 
relationship  may  be  listed  as  follows: 


v'y "vy  - t ■"  '/•>•  .-••Av'y- . - 


x^  p  cosp,  X£  ~  P  S  1  D(p j  X^  1 

“l  =  *1  cos91„  -  '3  s1n  ‘in  ; 
u 2  ""  *2  ’ 

u3  =  X1  Sin0in  +  x3  C0S8in  ; 

2  2  2 
%  =  ui  +  u2  ; 

e  =  tan  ^  (uju. )  ;  and 

u  '21'’ 

z  =  u~  . 
u  3 

Further,  an  expression  similar  to  Eq.  (9)  can  be  obtained  as 

2  2  2  2 
p  =  {p  cose,  cos*  -  z  sin  9.  )  +  p  sin  0 

u  '  in  in' 

$  =  tan  ^[(p  sino)/(p  cose.^  cos4>  -  z  sine.^)] 


z  =  p  sine.  cos<t>  +  z  cose. 
u  in  in 

Accordingly,  the  location  of  detector,  initially  identified  by  using  the 

parallel  coordinate  (p,e,z)  can  finally  be  rewritten  in  terms  of  the  unparallel 

system  (eu>ru,zu)  (see  Fig.  2~(8)  as, 

2  2  2.2 

p  =  (p.  cose.  cos<t>,  -  z.  sine.  )  +  p,  sin  <&,  , 

u'd  in  d  d  in'  d  d 

<j>  =  tan  lf(p.  sin4>.)/(p,  cose.  cos0,  -  z,  sine.  )]  and 

u  1 '  d  d  '  d  in  d  d  1  n  1 

z  =  p  ,  sine.  cosct>,  +  z.  cose.  .  (10) 

u  d  in  d  d  in 

Since  it  is  readily  found  that  the  geometric  configurations  of  the  inclined 
plume  vs.  the  coordinate  (y^y^y^)  equivalent  to  those  of  the  similar 
plume  with  axis  in  vs.  an  equally  inclined  detector  coordinate 
(UpU^uj),  some  pertinent  conversion  may  be  made.  Converting  the 
geometric  details  of  a  given  solid  angle  initially  represented  by  using  (9.0 
into  an  expression  in  terms  of  correspondi ng  unparallel  coordinate  (0U,S  ), 
the  results  obtained  for  the  case  of  mutually  parallel  coordinate,  Eq.  (9),  can 
directly  be  used.  For  this,  the  conversion  for  the  solid  angle  is  achieved  by 


11. 


8  =  COS  ^  r  -  S  i  ri  0  COSC  Sin  6.  +  CuStt  COS  6  .  1 

ij  1  in  i  n' 

C  =  tan  ^fsine  sinc/fsine  cose  cose.  +  cose  sine.  )1 
u  L  in  in'J 


(11) 


The  location  of  a  data  point  originally  identified  in  terms  of  the  coordinate  of 
plume  of  which  axis  is  inclined  by  an  angle  of  e.  with  respect  to  a 
perpendicular  coordinate  can,  therefore,  be  expressed  in  terms  of  a  spherical 
coordinate  centered  at  a  detector  on  the  combustor  wall  by  introducing  Eqs.  (9) 
and  (10)  i nto  Eq.  (11). 


( 2 )  Pi str i but i on-to-Di stribution  T ransformati ons 

WITHOUT  SWIRL  MOTION  -  As  explained  earlier,  the  details  of  i n-combustor 
spray  plume  may  be  obtained  in  an  equational  form,  e.g.  Eq.  (8),  expressed  with 
respect  to  the  nozzle  hole.  The  transformation  of  such  distribution  to 
individual  distributions  along  an  individual  optical  path  at  a  detector  is 
considered  for  a.  combustor  without  swirl  motion.  The  logical  step  for  this  goal 
would  have  been  introducing  the  correspond!' ng  terms  of  Eq.  (9)  into  Eq.  ;8). 

Such  transformation,  however,  results  in  an  extremely  complex  form  of  Eo.  (8) 
expressed  in  terms  of  the  spherical  coordinate  at  the  detector.  Furthermore, 
the  new  transformed  distribution  equation  becomes  too  difficult  to  be  used  in 
Eq.  (5).  In  order  to  alleviate  the  problem,  an  approximation  technique  was 
sought  in  such  a  way  as  to  conveniently  use  its  results  in  obtaining  the 
solution  of  Eq.  (5).  It  was  found  that  use  of  the  Taylor  expansion  method  in 
the  transformati on  not  only  greatly  facilities  the  steps  toward  the  final  goal 
but  also  lead  to  very  accurate  results. 
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p-terms : 

Letting  \  =  r  sine/p^,  A  5  cos(od  +  c, ) ,  with  f  s  p/p  one  finds 

7  1/2 

f(x)  =  (1  +  2AX  +  x  )  , 

f ' (x)  =  (A+x)/f  ,  and 
f"(x)  =  [f2(x)  -  ( A+x ) 2 ]/f 3 


For  \  =  x^  +  e  and  e<<1  the  above  may  be  rewritten  as, 

f(x)  =  f(Xj)  +  f’U-j)  e  +  f  •>{xi  )e2/2 

Further,  since  e  =  x_x^»  f(x)  rcay  be  expressed  as 

f(x)  =  p/pd  =  f0  +  fi  +  f2  x2  (12) 

where,  fQ  =  f(x1)  -  f’(Xj)  xi  +  f’Uj)  x-/2, 
f ' (xi )  -  f"(xi )  xi ,  and 

f2  =  nx^/2. 

For  x»l,  f(x)  =  p/pd  =  x  +  A 

$  terms: 

From  Eq.  (9),  having  f(x)  =  tan(o) 

=  (sin((>d  -  x  sine )/(coso>d  +  \  cose). 
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one  obtains 

f'(x)  =  -s i n ( <t>d  +  ?)/(cosod  +  X  cosOZ, 
f"(x)  =  -2  f '  ( x )  cosc/(cos<?d  +  x  cose), 

4>  ( x )  =  tan-1  ( f ) 
t>‘  (x)  =  f‘/(l+f2),  and 
♦"(x)  =  [-2  ff ‘ 2  +  (1+f2)  f"]/(l+f2)2. 

2 

Note,  for  cos$d  +  x  cose  =  0,  <t>(x)  =  */2,  <1* '  ( x )  =  -cos  e/sin(od+e)  and 
o" ( x )  =  0.  For  x  =  x.j+e  and  e<<l,  <p(x)  may  be  written  as  <p(x)  =  o ( x i )  + 
<t> '  (xi  )  £  +  )  e2/2. 

2 

In  view  of  Eq.  (8),  one  defines  4>(x)  =  <p  (x)  to  rewrite 
$(x)  =  «(x1)  *'  (x^  £  +  *"(x1)  £2/2 


where  <J>'  (x)  =  2<t(xi )  4>'  (x^ ) 

*"(x)  =  2<j,(Xi)  0“  ( x1 )  +  2  6 ' 2  ( x  n- ) 

Further,  since  e  =  x_X.j,  ^  may  be  expressed 

*(x)  -  <t,0  +  <5>i  X  +  $2  X2 

2 

where,  =  <f> ( x  1- )  -  (xi )  xi  +  x i /2 

^  =  'J* ( x i )  -  ,&"(xi )  xi 
*2  =  *"(xi)/2 

~  2 

For  x»l,  <p( x )  =  -e,  or  <*> ( x )  =  e  • 


(13) 


z-terms : 

Since  z  is  a  single  first  order  polynomial  function  of  r,  the  following  is 
obtained  from  Eq.  ( 9) . 
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z 


2 


2  2  2 
z  .  -  2  z  .  r  cose  +  r  cos  e 
a  a 


(14) 


Finally,  by  substituting  terms  in  Eq.  (8)  with  Eqs.  (12),  (13)  and  (14),  one 
obtains 


C  £  -[(r’-r  )/r  ] 

F  =  f  ,  e  L '  o'  wJ 
d 


where,  f  .  =  f  e  ^al  ^1  '^l  ^ , 


r  = 


b1/2  Cr 


r„  =  (ci 


■1/2 


a,  = 


b,  = 


c,  = 


a  pQ  +  b  $0  +  c  2/, 

-[a  +  b  4>^]  sine/pd  +  2  Czd  cose,  and 

2  2  2 

[a  p^  +  b  $2]  sin  0/pd  +  C-cos  e. 


(15) 


When  the  above  results  are  used  in  solving  Eq.  (5),  it  is  found  that  the  above 
approximation  method  generates  better  accuracy  when  the  integral  interval  is 
divided  into  smaller  subranges. 


WITH  SWIRL  MOTION  -  The  present  coordinate  transformati on  method  is  also 
considered  for  a  combustor  with  swirl  motion.  With  swirl  motion,  prior  to  the 
transformation,  the  spray  plume  distribution  equation,  Eq.  (8),  needs 
adjustments  to  take  into  account  the  bending  of  spray  plume.  Although  the  exact 
formation  mechanism  of  bent  sprays  in  the  presence  of  gas  swirl  in  the 
combustor  is  not  properly  understood  at  present,  the  following  two  probable 


cases  are  considered. 


Droplet  with  C onstant  Axial  Velocity : 

It  is  assumed  that  (1)  fuel  droplets  ejected  from  the  nozzle  hole  to  move  at 


a  constant  axial  velocity,  v^;  (2)  the  angular  velocity  of  the  swirl,  ft,  is 
constant  along  the  axis  of  the  plume;  (3)  the  angular  velocity  of  fuel  droplets 
along  the  plume  axis  is  the  same  as  Q.  With  the  above,  a  point  character i sti c 
at  a  location  ( p , <j> , Z )  of  no-swirl  condition  is  found  at  a  location  (p,$+<>  ,Z), 
due  to  the  swirl  during  a  period  of  time  At,  according  to 
4>s  =  (Q/vf)p, 

i 

«>s  =  ft/v^  and 

II 

<j>  =0. 

s 

2  2 

In  this  case  Eq.  (8)  becomes  F  =  Fq  exp[-ap  -  b(<p-<t>s)  -Cz  ].  The 
transformati on  of  the  above  equation  into  the  detector  coordinate  expression 
is  exactly  identical  except  for  a  new  set  of  equations  in  (j>-terms.  They  are: 

*(xi)  =  't'2  ( x  i ) , 

#'  (xi )  =  2->(xi )  <!>'  (xi ),  and 

♦"(Xi)  =  2-)>(xi)r(xi)'  +  24>,2(xi ) ,  (16) 


where,  *(Xi)  =  ♦(x1)  "  ♦$(xi). 

i 

<t»'  (x.j )  =  f '  (xi )  -  <t>s ( x i )  and 
)  -  4>s(xi ) . 


Droplets  with  Specified  Axial  Velocity: 

Melton  [21]  analyzed  a  case  for  the  plume  formation  as 


v 


f 


Vc/(  1  +  ctp ) 


where, v$  is  the  exit  velocity  of  the  fuel  droplet  at  the  nozzle  hole,  a  is 
about  0.085/Rs  and  Rs  is  the  radius  of  the  nozzle  hole.  Assumptions  of  the 
previous  case  are  also  similarly  taken  here  for  finding  <|>  =  ft-t,  since 

v ^  =dp/dt  =  vs/(l+ap)  and  p  +  ap^/2  =  vs  t. 

Combining  the  above,  one  finds 

<j>s  =  (ftp/vs)  (l+ap/2). 

Letting  <|)  ^  =  (Q/v  )  (l+a/2),  the  above  equation  may  be  expressed  in 

terms  of  0  n  to  obtain 
si 

<i>s  =  <t>sl  p(l+ap/2)/(l+o/2) , 

*)»$  =  <t>$^(l+«p)/( l+a/2 )  and 
<t>s  =  a/(l+a/2) . 

Substituting  respective  terms  in  Eq.  (16)  by  the  above  results,  the  coordinate 
transformation  of  bending  spray  plumes,  as  specified  here,  can  be  achieved. 

Regarding  the  plumes  with  their  axis  having  an  inclined  angle  of  9^  ,  the 
distribution-to-distribution  transformation  can  be  obtained  by  introducing 
Eqs.  (11)  and  (12)  into  Eq.  (15). 

COMPUTATION  OF  RADIATION  IN  DIESEL  COMBUSTION 

The  present  method  may  readily  be  used  for  analysis  of  radiation  heat 
transfer  processes  in  combustors  having  jet  flames.  An  analysis  may  be  made 
over  the  entire  wall  of  the  reactor  to  introduce  the  results  to  implement  a 
global  analysis  of  the  system  heat  transfer  including  convective  and  conduction 
heat  transfers,  which  will  enable  prediction  of  the  temperature  distribution 


17. 


into  and  over  the  wall.  When  an  unsteady  combustion  reactor  is  considered,  e.g. 
a  diesel  engine,  the  analysis  will  have  to  be  made  over  the  entire  cycle  period 
along  with  energy  conversion  process,  as  presently  attempted  in  our  laboratory. 
The  present  section,  however,  considers  the  computation  of  radiation  heat 
transfer  at  a  particular  engine  crank  angle.  Since  more  extensive  results  are 
reported  elsewhere  [22],  only  a  few  computations  are  described  here.  Among  the 
assumptions  employed  due  to  the  absence  of  either  sufficient  engine  data  or 
suitable  engine  modeling  are  (1)  the  values  of  constants  in  Eq.  (8)  are  for 
soot,  H^O,  CO^  and  burned  fuel/air  ratio,  8x10  0.01,  0.01,  and  1.0, 

respectively;  the  distribution  constants,  a,  b  and  c  are  0.6,  2.36,  and  2  for 
the  above  cases  accordingly.  {For  example,  the  CO^  distribution  in  the  plume 
may  be  written  as  F  =  {0.01  Exp  -  [0.6(p/R)  +  2. 43(0/^)  +  2.0(z/R)2]},  where 
R  is  the  piston  radius  and  is  ir/( number  of  plumes),  (2)  other  engine 
details  are,  the  surface  temperature,  500°K;  the  number  of  plumes,  4,  the  piston 
radius,  4.92  cm;  the  combustion  chamber  bowl  radius,  3.44  cm,  and  the  fuel 
composition  (3)  the  computation  of  the  temperature  distribution  is 

computed  by  using  the  equilibrium  adiabatic  temperature  relationship  to  the 
fuel/air  ratio  (e.g.,  using  a  new  method  [23]),  etc.  The  computed  radiation 
heat  flux  is  presented  presented  in  terms  of  the  following  nondimens i onal i zed 
terms:  the  normalized  spectral  energy,  Q  /Q  ,  and  the  normalized  directional 

A.  D  A. 

radiation  intensity,  I  ( p  ,<D,z)/Ib,  where,  Qbr  =  aT^,  Tr  =  2,400°K,  Ib  =  Q  /*, 
and  a  is  the  Stefan-Bol tzman  constant  (5.670xl0~®W/m2  K4). 

Shown  here  are  the  following  case  analyses:  (1)  the  spectral  radiation  heat 
incident  upon  locations  on  the  cylinder  head  along  the  axis  of  a  spray  plume  at 
r  =  0.0R,  0.5R  and  0.9R  (Fig.  3),  (2)  the  normalized  directional  radiation 
intensity  streaming  into  the  injection  nozzle  tip  (p  =  0.0,  0  =  0.0)  along  the 
direction  of  zenith  angle  in  the  hemispherical  volume  faced  by  the  tip  'Fig.  4) 


18. 


and  (3)  the  normalized  directional  radiation  intensity  along  the  locations  on 

the  cylinder  head  right  above  the  axis  of  the  spray  plume  in  directions  of 

zenith  angle,  9  =  17.6,  60.0  and  76.7  in  the  hemispherical  volume  faced  by  each 

location  (Fig.  5).  Without  discussing  the  details  of  the  results,  one  may  readily  find 

from  the  above  illustration  that  improved  radiation  heat  transfer  computations 

can  be  achieved  by  employing  the  present  new  coordinate  transformation  method 

in  its  modeling. 
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NOMENCLATURE 

A  band  absorptance 

a  constant  determining  species  distribution 
B  radiosity 

b  constant  determining  species  distribution 

c  speed  of  light  or  constant  determining  species  distribution 

e  Plank's  emission  power 

F  species  distribution  function 

f  soot  volume  fraction  or  transformation  function 

h  Plank's  constant 

I  radiation  intensity 

k  Boltzman  constant  or  imaginary  part  of  reflective  index 
n  real  part  of  the  reflective  index 
Q  black  body  radiation 
q  radiation  heat  flux 


R  radius  of  cyl i nder 
T  temperature 
t  time 


v  velocity 

x  component  of  cartesian  coordinate  X 
y  component  of  cartesian  coordinate  Y 

a  constant  for  plume  injection  speed 
S  reflectivity  of  cylinder  wall 
r,e,c  components  in  spherical  coordinate,  C 
e  emittance 

e  angle  of  optical  path  with  respect  to  normal 
k  volume  absorptance 
X  wavelength 

p,<i>,z  components  in  cylindrical  coordinate,  S 
t  extinction  coefficient 
$  transformation  function 
ft  angular  speed 

Subscripts 

b  blackbody 
d  detector 
g  gas 

in  inclined 
o  constant 
s  soot 


t  total 


2 


1 


u  unparallel 
\  spectral 

es  optical  path  of  solid  angle 
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Detector  location  with  respect  to  flame  plume. 

(A)  Parallel  Coordinate  and  (B)  Unparallel  Coordinate 
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Directional  radiation  intensity  at  the  injection  nozzle  streaming  through 
individual  direction  of  zenith  angle. 
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Flow  Chart  of  HEATRAN 
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C***  HEAT  TRANSFER  MODEL  OF  DIESEL  COMBUSTION 
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IMPLICIT  REAL*8  (A— H ,0—2) 

COMMON  /BLK/RPM,RP,RL,RCS,CR,CM,HN,QC 

COMMON  /BLKR1/PI ,PIST,PlL 

COMMON  /BLKR2/CFRE (hO) , ERF C , ERFCL 

COMMON  /BLKR6/S0CF  (3) ,SOCG  (3)  .SOCWN  (3)  , VLGT 

COMMON  /BLKR7/REF L , REFTP, REFQ.REFW, REF WL, CVLWL 

COMMON  /BLKR8/0L,PfPE,T,T0,W 

COMMON  /BLKR9/AE (20) 

COMMON  /BLKV2/R  (20)  , QV  (20) 

COMMON  /BLKT6/TM(360) 

COMMON  /BLKC1/XC  (60)  ,YC  (60)  ,C  (60)  ,CK  (60) 

COMMON  /BLKC2/A (60,60)  ,B  (60,60) 

COMMON  /BLKC3/D  (60)  , UBD  (60)  , UO  (60) , UB  (60) 

COMMON  /BLKCVREFT.REFOT.REFC.REFK.REFLC.NPT.NU.MBC.NBCT, 
COMMON  /BLKC5/L1  (80)  .  L2  (80)  ,  L3  (80)  .  I  1  (80)  .  I  2  (80)  ,  I  3  (80)  , 

1  Ji  (8o)  ,J2  (80)  ,J3(80) 

DIMENSION  UB1  (60)  , QTT  (11) 


INPUT  AND  OUTPUT  FORMATS 


FORMAT (2F7.1 .3E10.3.5I 5) 
FORMAT  (2F10.4,2E10-3,F10.U) 
FORMAT  (915) 

FORMAT  (FI  0.5) 

FORMAT  (F9. 3, F12.4) 

FORMAT  (FI  0.1) 

FORMAT  (7F7-1) 

FORMAT  (3 (8E9 . 3 . /)  .4E9-3) 

INPUT  COMMON  PARAMETER 


READ  (5.33) RPM,RP,RL,RCS,CR,CM,HN,QC 
FORMAT (F7.1.3F9.3.3E9.1.F11.1) 


CONSTANTS 


P I “4 . D0*DAT  AN (1 .00) 
PIST-DSQRT(PI) 
p  |  i,.p  |  *p  |  *p  |  *p  | 

CALL  I  CFRE  (40) 


R£FQ=1 .880502 
REFTP=2.4D3 
REFW=1 .25015 
REFWL=6 . DO 

HCK=REFWL*REFTP*1  .D-4 

Disperion  constants  of  soot  optical  constant  from  Lee  £  Tien. 

SOCF  (1) -7.2026 
SOCF  (2)=4. 07027 
SOCF  (3) =4.47D28 
SOCG  (1) =1 .2015 
SOCG  (2) =5-901 5 
SOCG  (3) =5.6015 
SOCWN (1)  =.0D0 
SOCWN  (2) =1 .25015 
SOCWN  (3) =7-25015 
DO  100  1-1,3 

SOCF  ( I ) =S0CF  (I) *3.183D3/R£FW/REFW 
SOCG  (l)-SOCG  (D/REFW 
SOCWN ( I ) =S0CWN ( I )  /REFW 
VLGT=2. 998401 4/REFWL/REFW 
T0-1 .D2/REFTP 
OMEG=.53D2 
TS=.5D3 

Finite  Difference  constants  of  Conduction 


READ  (5.1) REFT, REFDT.REFC. REFK, REF LC.NPT, NLP, NU, NBC, NBCT 
READ  (5 ,2)  (XC  (N)  ,  YC  (N)  ,  C  (N)  ,  CK  (N)  ,  UO  (N)  ,  N-1  ,  NPT) 

READ  (5,3)  (LI  (I)  ,11  (I)  ,J1  (I)  ,L2(I)  ,12(1)  ,J2(I)  ,L3(I)  ,I3  0: 
+  J3(l)  ,  1=1, NLP) 

C  A  A  A  A  A  A  A  A  A  it  it  A  A  A  it  A  A  A  A  A A  A  *  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A  A A  A  A  A  A  A  A  A  A  A  it  A  A  A  A  A  A  A  A  A : 
C  AAA 

C**A  MA I N  PROGRAM 


■A  AAA  AAA 


CALL  MATRIC 
CALL  THERMO 
DO  200  K=1 0,50,10 
TG=TM  (K) 

CALL  CONV (OMEG ,TG ,TS) 

WRITE (6,5)  ( (R ( I )  , QV ( I ) ) ,1=1,11) 
WRITE  (6,6)  TM  (K) 

CALL  RANDO(CAO) 

WRITE  (6,7) CAD,  (AE  (I ) ,1=1,11) 

DO  120  M=1 ,11 

QTT  (M) =QV (M) /I . 0D3+AE (M) *1 88 . 05DO 
IF  (K.LE.10)  GOTO  160 
READ (5,8)  (UB1  (L)  , L  =  1 , NBC) 

DO  140  1=1,11 
UB1  ( I +1 4)  =QTT ( I ) 

TK0= (K-2)  *1 .DO/3-OD2 
TK= (K-1)  *1 .00/3.002 
CALL  COND (UB1 .TKO.TK) 


SUBROUTINE  THERMO 
IMPLICIT  REAL*8  (A-H ,0-2) 

COMMON  /BLKT3/PG  (360)  ,WG  (36O) 

COMMON  /BLKT4/PA , CA1 .DCA.TA I R.RHM 
COMMON  /BLK/RPM, RP , RL , RCS , CR , CM, HN , QC 
.  COMMON  /BLKT5/XA I R  (4)  , XR  (4) , XO  (4)  , ZNA I R , ZNO 
-  COMMON  /BLKT6/TM(360) 

COMMON  /BLKT7/VMAXS  (360) 

COMMON  /BLKT8/SKG  (360,2) 

COMMON  /BLKT9/E  (360) 

D I  MENS  I  ON  CAT  (36O)  , QX  (360)  , DQX  (36O)  , XC02T  (36O) 

DIMENSION  PS  (360,2)  , QS (360,2) , DQS  (360,2) 

DIMENSION  PERS  (2,5,5)  ,ZN  (4)  ,CV(4)  , U  (4)  ,  I  DATE  (3.5)  ,  I  CYCLE  (5) 

1  FORMAT(6F9.3,F11  .1) 

2  FORMAT  (I  4) 

3  FORMAT  (5F9 - 3. 14,312) 

4  FORMAT  (10F8. 2) 

11  FORMAT  (T5, ' Eng i ne  parameter  :',/ ,T8 ,’ crank  arm  -'.T28.F6.3, 

1*  cm’ ,/,T8, 'connecting  rod  - ' ,T28, F6 . 3, '  cm ' , / ,T8 , 1 cy 1 i nder  ' 
2’diemeter  - ' ,T28 , F6 . 3 • 1  cm ',/ ,T8 , 1  compress i on  ratio  -',T28, 
3F6.2,//,T5, 1  Fuel  properties  :',/,T8,'C:H  -  1 ,T28, I  2 , 1 : 1 , I  2 , / , 
4T8,'low  heat  value  -  1 ,T26 , F9 . 0 , '  KJ/Kmole',/) 

12  FORMAT ('1 ', ’TEST  : 1 , I  4 ,T20 , ' Date  :  ' , I  2. '/' , I  2. '/'  .  I  2.//.T5, 

I  1  inlet  press.  .T30.F6.1 ,  '  Kp' ,//,T5,  '  inlet  temp.  -\T30, 
2F6.1,'  K' ,//,T5, 'relative  humidity  - ' ,T30, F6 . 1 , '  %',//,T5, 

3 1  initial  sampling  - ' ,T30 , F6 . 1 , '  deg .', //,T5 ,' samp  1 i ng  ', 
4'interval  -',  T30.F6.3,'  deg  . ' , //,T5 , ' eng i ne  speed  -', 

5T30, F6.0,  1  rpm'  ) 

F0RMAT(/,'  CYCLE  : ' , I 4.T30, ‘ there  are',14,'  pressure  data  ', 

I I  in  the  cycle. ' ,T80, ' UNIT  :  CA=deg . ; P=Kp ' , // ,T4 , ' CA ' , T71 . 
2T19, ' CA 1 ,T26, ' P ' ,T34, 'CA' ,T41, 'P' ,T49, 'CA‘ ,T56, ' P ' ,T64, ■ CA ' , 
3T71  , 1 P’ ,T79, ' CA ' ,T86, ' P ' ,T94, 1 CA ' ,T101 , ' P ' ,T109, ' CA ' ,T1l6, 
4'P',T11,'P') 
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FORMAT (8(2X,F5.1,2X,F6.1)) 

FORMAT (F5.0.6F11 -3) 

FORMAT (F5.0.3F11 .3) 

FORMAT  (/ ,T5 , '  mean  effective  press.  - ' ,T30, F6 . 1 , 1  Kp1,/, 

1T5,'ind.  horse  power  -  1 ,T30 , F6 . 2 , '  HP  1 , /,T5, 1 f ue 1  air  ratio', 

2'  ,T30,F7. 5,/. T5, ' thermal  efficiency  - 1 ,T30 , F6 . 1 , '  %',/, 

3T5,'deg.  of  const,  vol.  -  1 ,T30, F6. 1 , 1  %',/) 

FORMAT  (/,T5, 'max .  pressure  1 , F  7 . 1 . 1  Kp  at  1 . F  5  - 1 , '  deg.  CA',/, 
1T3,'max.  heat  release  rate  1  ,  F  7  -  ^ .  1  kw/deg.  at  ' , F  5  - 1 . 

2' deg.  CA  * . /) 

FORMAT  ('1 ' ,T20, 'NORMALIZED  PRESSURE  vs  CRANK  ANGLE ', /.T20 , 
Preference  pressure  -  p"  (ps  i )  1  ,/,T6 ,  ‘  P/P"  1  ,T1  7 ,  '  CA  1 ) 

FORMAT (' 1 1 ,T20, 'NORMALIZED  RELEASED  HEAT  vs  CRANK  ANGLE1, /,T20 
1, 'reference  released  heat  Q"  (KJ)  1  , /,T6 , 1  (J/Q"  '  ,T1  7 ,  '  CA ' ) 

FORMAT  (' 1‘ ,T20, 'NORMALIZED  HEAT  RELEASED  RATE  vs  CRANK  ANGLE',/, 
1T20 ,  1  ref erence  heat  release  rate  -  DQ" (KJ/deg . ) ' , / ,T4 , ' DQ/DQ" ' ) 
FORMAT  (' 1 1 ,T7 , 1  crank  ang 1 e ' ,T25 pressure ' ,T46 ,' heat  release', 
1T60, 1  heat  release  rate ', /,T1 0 , 1  deg . 1 ,T27 ,' Kp ' ,T50, ' KJ 1 ,T65 , 

2 1 KJ/deg .CA 1 , / ,T22 . 1 ave . 1 ,T32 , 1 S . D . ' ,T42 , ' ave . 1 ,T52 , ' S . D . ' , / , 
3T64, 'ave. 1 ,T74, 'S.D. ') 

FORMAT  (T10.F5.1  ,T20,F7.1  .T30.F6.1  (4X.F7.4) ) 

FORMAT  (' 1 1 ,/,T3, 'The  engine  performance  in  the  tests  :',/, 

IT!*,  'Date'  ,T13,  '  Cyc  1  es  1  ,T29,  '  imep '  ,T46,  1  i  nd .  HP  '  ,T69,  1  F/A  '  , 

2T86 , ' ef f i c i ency ' , T1 03 , ' deg .  of  const,  vol . 1 ,/,T25. ' ave. 1 ,T35. 
3'S.D. ' ,T45, 'ave. ' ,T55. 'S.D. ' ,T65, 'ave. 1 ,T75. ’S.D. ' ,T85, 'ave. ' , 
4T95. 'S.D. ' ,T105, 'ave. ' ,T115, 'S.D. ' ,/) 

FORMAT (2X, I  2 , 1 / ' , I  2 , 1 / 1 , I  2 , 4X . I  4 , 3X , 5 (F 1 0 . 2 , F 1 0 . 3) ) 

MC=CM 

NH=HN  ;  . 

WRITE (6,11)  RCS , RL , RP , CR , MC , HN , QC 

ITEST=0 

READ  (5.2)  INDEX1 
I NDEX1= I NDEX1 -1 
I F  (  INDEX1  .LT.  0  )  GOTO  800 
I TEST= I TEST+1 
I  CYCLE (ITEST) =0 

READ (5,3) PA ,TA I R , RHM, CA1 ,DCA, I  SAMP,  (I  DATE  (M, ITEST)  ,M=1 ,3) 

WRITE  (6,12)  ITEST,  (I  DATE  (M, ITEST)  ,M=1 ,3)  , PA ,TA I R , RHM, CA1 , DCA , RPM 

READ  (5,2)  INDEX2 

I NDEX2= I NDEX2-1 

I F  (  INDEX2  .LT.  0  )  GOTO  700 

DO  350  1=1  , I  SAMP 

E  ( I ) =CA1+  ( I  - 1 ) *DCA 

I  CYCLE  (I TEST) = I  CYCLE  (I TEST) +1 

READ  (5,4)  (PG  ( I )  ,1=1,  I  SAMP) 

WRITE  (6,13)  ICYCLE  (ITEST)  , I  SAMP 
WRITE  (6,14)  (E  (J)  , PG  (J)  ,J  =  1 ,  I  SAMP) 

CALL  PERFOR ( I  SAMP , P I  ME , HP  I , F ATH , EF F , DCV) 

CAT  (1)  =CA1 
DO  400  J=2 , I  SAMP 

CAT  (J) =CA1  +  (J-1) *DCA 
CALL  VOLUME  (CAT  (J)  ,V,DV) 

TM  (J) =PG  (J) *V/  ( (ZN0+ . 25D0*HN*E  (J-1 ) ) *8 . 31 4D6) 

CALL  SPHEAT  (TM  (J)  ,CV,U) 

E  (I SAMP+1) =0.D0 


OQX  (J)  =  {E  (J+1)  -E  (J-1) )  *QC/2  . DO*DCA 

QX  (J)  =E  (J)  *QC 

ZNT=ZNO* (1  .+  .250O*HN*E (J) ) 

XC02T  (J)  =  (ZNO*XO  (3)  +HN*E  (J)  )  /ZNT 
400  WHITE  (6,15)  CAT  (J)  ,  PG  (J)  ,V,TM(J)  ,  DQX  (J)  ,QX  (J)  ,  XC02T  (J) 
DO  500  J=51,210 

CAT  (J)  =  (CA1+ISAMP*DCA)  +  (J-51) *DCA 
ZNL=ZNO*  (1 .D0+.25D0*HN*E (I  SAMP)  ) 

ZN  (1 ) =ZN0*X0  (1) -  (CM+.25D0*HN) *E  (I  SAMP) 

ZN  (2) =ZN0*X0  (2) 

ZN (3) =ZNO*XO  (3) +CM*E  ( I  SAMP) 

ZN  (4)  =ZNO*XO  (4)+.5D5*HN*E  (I  SAMP) 

CALL  SPHEAT  (TM(J-I)  ,CV,U) 

CALL  VOLUME  (CAT  (J)  ,V,DVL) 

ACV=CV  (1)  *ZN  (1)  +CV  (2)  *ZN  (2)+CV  (3)  *ZN  (3) +CV  (4)  *ZN  (4) 
TM  (J) =TM  (J-1 ) -PG  (J-1) *DVL*DCA/ACV/1 .D6 
PG  (J) =ZNL*8 . 31 4D6*TM  (J) /V 
500  WRITE  (6,16)  CAT  (J)  , PG  (J)  ,V,TM(J) 

WRITE  (6,17) PI  ME, HP  I .FATH.EFF.DCV 
600  I CT=I CYCLE  (ITEST) 

CALL  STATIS (PERS  (1 , 1 , I CT) , PERS (2 , 1 , I CT) , I CT, P I  ME) 

CALL  STATIS (PERS (1 , 2 , I CT) , PERS (2 , 2 , I CT)  , I CT, HP  I ) 

CALL  STATIS  (PERS  (1,3,1 CT)  ,PERS  (2,3, I CT) , ICT.FATH) 

CALL  STATIS  (PERS  (1 ,4,  I  CT)  .P'ERS  (2 , 4 ,  I  CT)  ,  I  CT ,  E  F  F) 

CALL  STATIS  (PERS  (1 ,5. I CT)  .PERS (2,5, I CT) , I CT.DCV) 

DO  620  L*1 , I  SAMP 

CALL  STATIS  (PS  (L , 1 )  , PS  (L , 2)  , I CT , PG  (L) ) 

E(L)=E(L)*QC 

620  CALL  STATIS(QS(L,1)  ,  QS  (L  ,  2)  ,  I  CT,  E  (L) ) 

DO  630  1=1 , I  SAMP 
63c  VMAXS ( I )  =PG ( I ) 

CALL  VMAX ( I  SAMP , PMAX , LPMAX) 

CAPMA  l=CA1+  (LPMAX-1 ) *DCA 

DQS  (1 , 1)  *0 .  DO 

DQS  (1 , 2)  -0 .  DO 

DQS  ( I  SAMP , 1 ) =0 . DO 

DQS  ( I  SAMP , 2) =0 . DO 

I SAMP1=I SAMP-1 

DO  640  L=2 , ISAMP1 

DQ= (E (L+1)  -E  (1-1) ) /  (2 . DO*DCA) 

I F (  DQ  . LT .  0 . DO  )  DQ-0.00 
WG  (L) =DQ 

640  CALL  STATIS  (DQS  (L , 1)  , DQS  (L, 2) , ICT.DQ) 

WG  (1 )  *0 .  DO 
WG  (ISAMP) =0 .DO 
DO  650  1=1,1  SAMP 
650  VMAXS ( I ) =WG ( I ) 

CALL  VMAX  ( I  SAMP , OQMAX , LDQMAX) 

CADQM=CA1+  (LDQMAX-1) *DCA 
WRITE  (6,18) PMAX , CAPMAX , DQMAX , CADQM 
GO  TO  300 
700  WRITE (6,19) 

DO  710  1*1 , ISAMP 
SKG  (I  ,1)=PS  (I  ,1) 

710  SKG  (1,2)  =PS  (1,2) 
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CALL  SKETCH  (I  SAMP, DCA.CA1) 

WRITE  (6,20) 

DO  720  1=1 , ISAMP 
SKG  (1,1) =QS  (1,1) 

720  .  SKG  (I  ,2)=QS  (I  ,2) 

CALL  SKETCH (I  SAMP, DCA.CA1) 

WRITE  (6,21) 

DO  730  1=1 , ISAMP 
SKG (I , 1 ) =DQS  (I ,1) 

730  SKG  (1,2)  =DQS  (1,2) 

CALL  SKETCH (I  SAMP, DCA.CA1) 

WRITE  (6,22) 

DO  750  1=1 , ISAMP 
CA=CA1+  ( I  —  1 ) *  D  C  A 

WRITE  (6,23)  CA.PS  (I  ,1)  .PS  (I  ,2)  ,QS  (I  ,1)  ,QS  (I  ,2)  , 

+  DQS  (1,1)  ,DQS  (1,2) 

PG  (I)  -PS  (1,1) 

750  WG ( I ) -DQS  (1,1) 

DO  7&0  1=1 , ISAMP 
760  VMAXS  ( I ) =PG  ( I ) 

CALL  VMAX  ( I  SAMP , PMAX , LPMAX) 

DO  770  1=1 , ISAMP 
770  VMAXS  (I)  =WG  (I) 

CALL  VMAX  (ISAMP, DQMAX.LDQMAX) 

CAPMAX=CA1+  (LPMAX-1 ) *DCA 
CADQM=CA1+  (LDQMAX-1) *DCA 
WRITE  (6,18) PMAX , CAPMAX , DQMAX, CADQM 
GO  TO  100 
800  WRITE  (6,23) 

DO  850  K=1 , ITEST 

850  WRITE  (6,25)  (I  DATE  (M, K)  ,M=1 ,3)  , I  CYCLE (K) ,  ( (PERS ( I  , J , K) , I =1 , 2)  , 


U-1,5) 

RETURN 

END 

C»> - <« 

c»>  «< 

C»>  SPECIFIC  HEAT  AND  INTERNAL  ENERGY  «< 

C»>  .  «< 

C»> - «< 


SUBROUTINE  SPHEAT  (T, CV, U) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  C \l(k)  , U  (4) 

1  FORMAT  (T20, 1 ERROR:SPHEAT  -  Temp .  =  ' , El  2 .5) 
IF  (  T  .GT.  0 . DO  )  GOTO  200 
WRITE  (6,1) T 
STOP 

200  T18»T*1.8D0 

T18ST=DSQRT  (Tl8) 

Tl8SQ=Tl8*Tl8 
Z1 9=1 .986DO 
Z83=8 . 31 43DO/Z1 9 
Z8318=Z83/1 .800 
Z5b-5.kD2 
Z5^ST=DSQRT(Z514) 

Z5ASQ=Z54*Z54 
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Z1 8LG=DL0G  (T18/Z5M 

CV (1 )  =  (1 -1515D1-1 . 72D2/T1 8ST+1 -53D3/T1 8-Z1 9)  *Z83 

cv (2) ={9.4700-3.4703/118+1 . 1 606/T1 8SQ-Z1 9) *z83 
CV (3)  =  (1  .6201-6. 52D3/T18+1 .41D6/T18SQ-Z19) *Z83 
CV (4)  =  (1 .998D1-5.6702/T18ST+7.5D3/T18-Z19) *  Z  8  3 
u  (1)  =  (9.528DO*  (T18-Z54) -3.44D2* (T18ST-Z54ST)+1 .5303*21810  *Z83l8 
U (2)  =  (7-48300* (T18-Z54) -3 . 47D3*Z 1 8LG-1 .16D6/TI8+ 

1  1 .16D6/Z54) *Z8318 

U  (3)  =  (14.21300*  (T18-Z54) -3.47D3*Zi8lg-i . 4 1 D 6 * 

1  (1 .D0/T18-1 .D0/Z54) ) * Z 8  3 1 8 

U  (4)  =  (1.787301*  (T18-Z54) -1.1 94D3*  (T18ST-Z54ST) 

1  +7.5D3*Zl8LG)  *Z83l8 

IF  (  T18  .LT.  4.D3  )  GOTO  900 
CV  (1)  =  (CV  (D  +  .5D-5*  (T18-4 .03) )  * Z 8 3 

U  (1 )  =  (U  (0+.5D-4*  (.500*  (T 1 8SQ-Z54SQ)  -4.D3*  (Tl8-Z54) ) )  *Z83l8 
900  RETURN 
END 


C»> - «< 

c»>  «< 

C»>  VOLUME  OF  COMBUSTION  CHAMBER  «< 

C»>  «< 

c»> - <<< 


SUBROUTINE  VOLUME (CA.V.DV) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR1/P I , P I  ST, P I  4 

COMMON  /BLK/RPM , RP , RL , RCS , CR , CM, HN , QC 

100  I F  (  CA  .GE.  0 . DO  )  GOTO  200 
CA=CA+3.6D2 
GO  TO  100 

200  IF  (  CA  .LE.  3.6D2)  GOTO  400 
CA=CA-3.6B2 
GO  TO  200 

400  THETA=CA*PI/1  .8D2 
A=P I *RP*RP 

VC=A*2 . DO*RCS/  (CR-1 .DO) 

XI -DSQRT  (1 .DO-  (RCS*DS I N  (THETA) /RL) **2) 

V- (RCS* (1 .DO+DCOS (THETA) ) +RL* (1 .D0-X1) ) *A+VC 
DV= (RCS*RCS*DS I N  (2 . DO*THET A) /  (2 . D0*RL*RL*X1 ) -RCS* 
1  DS IN  (THETA)) *A*P 1/1 .8D2 
RETURN 
END 


C»> - «< 

C»>  «< 

C»>  PISTON  WORK  «< 

C»>  «< 

c»> - <<< 


SUBROUTINE  WORK  (CA1 ,DCA,NO) 
IMPLICIT  REAL*8  (A-H.O-Z) 
COMMON  /BLKT3/PG  (360) ,WG (3&0) 
WG  (1)»0.D0 
DO  100  J-2.N0 

CA-CA1+  ( J  —  1 ) *DCA 
CAO-CA-DCA 

CALL  VOLUME (CAO.VO.DVO) 
CALL  VOLUME (CA.V.DV) 
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DW=  (PG  (J)  *DV+PG  (J-1)  *DV0)  *DCA/2.D6 
100  WG (J) =WG  (J-1) +0W 
RETURN 
END 


C»> - <« 

c»>  «< 

C»>  AIR  COMPOSITION  «< 

C>>>  <<< 

c»> - <<< 


SUBROUTINE  A  I R (TA I R , PA , RHM) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKT5/XAIR  (L)  , XR  (4)  , XO  (4)  .ZNAIR.ZNO 
DIMENSION  B(5)  ,F  (5) 

DATA  F/.17803D0, .2563D0, .363200, .507300, .6988DO/ 
R=(1 .8D0*TAIR-5.3D2)/1 .01 
B (1)  =R* (R*R— 1 .DO) *  (R-2.00) /2.UQ1 
B  (2)  =  -R*  (R-1  .DO)  ft  (RftR-4  .  DO)  /6.D0 
B  (3)  *  (RftR-1  .DO)  *  (RftR-4  .DO)  /h  .DO 
B  (1*)=-R*  (R+1  .DO)  ft  (RftR-L  . DO)  /6.DO 
B  (5)  =R*  (RftR-1  .DO)  *  (R+2 .00)  /2  . i*D1 
PVS-O.DO 
DO  100  1=1,5 
100  PVS-PVS+B (l)  *F (I) 

PV=RHMftPVS/1 .02 
ASH=PV/ (PA-PV) 

XA I R  (1 )  =  . 2099D0/ (1 .DO+ASH) 

XAIR  (2) =.7898D0/ (1 .DO+ASH) 

XA I R (3) = . 30-3/ (1 -DO+ASH) 

XAIR  (A) -ASH/  (1 . DO+ASH) 

RETURN 
END  . 


C»> - «< 

C»>  «< 

C>»  FUEL  CONSUMPTION  <« 

C»>  «< 

c»> - <<< 


SUBROUTINE  EXTENT  (TR ,  I S AMP , TBDC , T 1 ) 

IMPLICIT  REALft8  (A-H.O-Z) 

COMMON  /BLKT3/PG  (360)  ,WG (360) 

COMMON  /BLKT9/E  (360) 

COMMON  /BLKTL/PA.CA1  .DCA.TAIR.RHM 
COMMON  /BLK/RPM,RP,RL,RCS,CR,CM,HN,QC 
COMMON  /BLKT5/XAIR  (1*)  ,XR(U)  ,X0(4)  .ZNAIR.ZNO 
DIMENSION  C V  (k)  ,U(U) 

BETA-1  .  OO/'CR 

TBOC-  (TAIR+BETAftTR) /  (1 . DO+BETA) 

CALL  VOLUME (0 . DO , VO . DVO) 

ZNO=PA*VO/  (8.31  A*D6*TBDC) 

ZNAIR-ZNO/ (1 .DO+BETA) 

ZNR-ZNA I R*BET  A 

CALL  VOLUME (CA1 ,V1 ,DV1) 

T1-PG  (1) ftVI/  (8.31 UD6*ZN0) 

CALL  SPHEAT  (T1 , CV, U) 

U1-0.D0 
DO  100  K-1 ,4 


8 


XO (K)  =  {ZNA I R*XA I R (K) +ZNR*XR (K) ) /ZNO 
100  U1=U1+ZN0*X0  (K) *U  (K) 

CALL  WORK (CA1 , DCA , ISAHP) 

E  (1 )  =0  .  DO 
DO  4.00  J  =  2 ,  ISAHP 
CA=CA1+(J-1) *DCA 
CALL  VOLUME (CA.V.OV) 

T=PG  (J) *V/  (  (ZNO+. 2500*HN*E (J-1) ) *8.31406) 

CALL  SPHEAT  (T, CV, U) 

E  (J)  =  (ZNO*  (XO  (1)  *U  (1)+X0  (2)  *U  (2)+X0  (3)  *U  (3)+X0  (4)  *U  (4) )  -U1  + 
1  WG (J) ) / (QC+ (CM+. 25DO*HN)  *U (1 )  -CM*U (3) - . 5D0*HN*U (4) ) 

IF  (  E (J)  .LT.  O.DO  )  E  (J) =0.D0 

I  F  (  E  (J)  .GT.  ZN0*X0(1)/(CM+.25D0*HN)  ) 

+  E  (J)  =ZNO*XO  (1)  /  (CM+.25D0*HN) 


400  CONTINUE 
RETURN 
END 

C»> - «< 

c»>  <« 

C»>  ENGINE  PERFORMANCE  «< 

C»>  «< 

C»> - <<< 


SUBROUTINE  PERFOR ( I S AMP , P I  ME , HP  I , F ATH , EFF , DCV) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  ZN  (4)  ,  U  (4)  ,  CV  (4) 

COMMON  /BLKT3/PG  (360)  ,WG  (360) 

COMMON  /BLKT9/E  (360) 

COMMON  /BLKT4/PA.CA1 , DCA , TA I R , RHM 
COMMON  /BLK/RPM . RP , RL , RCS , CR , CM , HN , QC 
COMMON  /BLKT5/XA I R  (4)  ,XR(4)  ,X0(4)  .ZNAIR.ZNO 
1  FORMAT  (/,T5, ' K=' , I 3,/,T5, 'The  assumed  residual  gas  temperatu', 
1 're  ',' TR= 1 , F6 . 1 ,/ ,T30 compos i t i ons  X02= 1 , F 6 . 4 , '  XN2=',F6.4, 

2 ' XCQ2 ' , ,F6.4, '  XH20= ' ,F6.4,/) 

100  K=0 

CALL  AIR  (TAIR.PA.RHM) 

200  TR=TA I R 

DO  250  1=1,4 
250  XR  ( I ) =XA 
300  K=K+1 

320  WRITE  (6,1) K,TR,XR 

350  CALL  EXTENT  (TR, I SAMP.TBDC.T1) 

400  ZNL=ZNO*  (1 .D0+.25D0*HN*E (I  SAMP)) 

ZN  (1 ) =ZNO*XO (1 ) - (CM+.25D0*HN)  *E (I  SAMP) 

ZN  (2) =ZNO*XO  (2) 

ZN (3) =ZNO*XO (3) +CM*E ( I  SAMP) 

ZN  (4)  =ZNO*XO  (4)+.5D0*HN*E  (I  SAMP) 

CAL*CA1+ (I  SAMP-1) *DCA 
CALL  VOLUME  (CAL.VL.DVL) 

TL=PG ( I  SAMP) *VL/  (8 . 31 406*ZNL) 

DCAL=  (3 • 6D2-CAL) /I .02 
CALL  SPHEAT  (TL.CV.U) 

UL=U  (1)  *ZN  (1 )  +U  (2)  *ZN  (2)  +U  (3)  *ZN  (3)+U  (4)  *ZN  (4) 

PL=PG  ( I  SAMP) 

DO  450  1=1,100 

CALL  SPHEAT  (TL.CV.U) 
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ACV=CV  (1 )  *ZN  (1 )  +CV  (2)  *ZN  (2)  +CV  (3)  *ZN  (3)  +CV  (4)  *ZN  (4) 

TL=TL-PL*DVL*DCAL/ ACV/1 .D6 

CA=CAL+I*DCAL 

CALL  VOLUME  (CA.VL.DVL) 

450  PL=ZNL*8.314D6*TL/VL 

I F  (DABS  (TR-TL) -5. DO)  600,600,460 
460  I F  (K-10)  500,500.600 
500  TR=TL 

DO  550  1=1,4 
530  XR  (I) =ZN  (I) /ZNL 
GO  TO  3C0 

600  IT=  ( . 1 8D3-CA1 ) /DCA 
CATDC=CA1+ I T*DCA 

WTDC=WG ( I T)  +  (WG  (IT+1) -WG (IT) ) *  ( . 1 8D3"CATDC) /DCA 
UB=ZN  (1 )  *U  (1 )  +ZN  (2)  *U  (2)  +ZN  (3)  *U  (3)  +ZN  (4)  *U  (4) 

CALL  VOLUME  (0 . DO , VB , DVB) 

CALL  VOLUME (.18D3.VT, DVT) 

P I ME= (WG ( I  SAMP) -WTDC+UL-UB) / (VB-VT) *1 . D6 

EFFV= (1 .DO-1 .DO/CR**.4DO) *1 .D2 

FATH= (1 .2D1*CM+HN) *E  (I  SAMP) / (2 . 8967DI *ZNO) 

CALL  SPHEAT  (TAIR.CV.U) 

UTBDC=XO  (1)  *U  (1)+X0  (2)  *U  (2)+X0  (3)  *U  (3)+X0  (4)  *U  (4) 

CALL  SPHEAT  (T1 ,CV,U) 

UT 1  =X0  (1)  *U  (1)+X0  (2)  *U  (2)  +X0  (3)  *U  (3)+X0  (4)  *U  (4) 
WALL=ZNO*  (UTBDC-UT 1 ) +WG ( I  SAMP) +UL-UB 
HP  I =WALL*RPM/  (.745700*6. D1) 

QALL=QC*E  (I  SAMP) 

EFF=WALL*1 . D2/QALL 
DCV=EFF*1 .D2/EFFV 
RETURN 
END 


C>» - <<< 

C»>  «< 

C»>  STATISTICAL  ANALYSIS  «< 

C»>  «< 

c>» - <<< 


SUBROUTINE  ST AT  I S  (AVE 1 . SD1 , N , SAMP) 

IMPLICIT  REAL*8  (A-H.O-Z) 

IF  (  N  .GT.  1  )  GOTO  100 
AVE 1 =S AMP 
SD1=0.D0 
GO  TO  900 

100  AVE2=  (N— 1 ) *AVE1 /N+SAMP/N 

SD1-SD1**2*  (N-2)  /  (N-1) 

SD1=SD1+(AVE1-AVE2) **2/ (N-1 )  +  ( (SAMP-AVE 2) **2)  / (N-1) 
I F  (SD1)  200,400,300 
200  SD1=DABS  (SOI) 

300  SD1=DSQRT (SD1) 

400  AVE1*AVE2 
900  RETURN 
END 


C>» - <<< 

c»>  «< 

C»>  THE  MAXIMUM  OF  A  SERIES  «< 

C»>  <« 


10 


c»> 


«< 


SUBROUTINE  VMAX (NO, SMAX, LSMAX) 
IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKT7/S  (360) 

LSMAX=1 
SMAX=S (1) 

IF  (  NO  .LE.  1  )  GOTO  900 
DO  100  1=2, NO 

IF  (  SMAX  .GE.  S(l)  )  GOTO  100 
SMAX=S  (I) 

LSMAX= I 

100  CONTINUE 
900  RETURN 
END 


C>» - <<< 

c»>  .  «< 

C»>  SKETCH  «< 

C»>  «< 

c»> - <<< 


SUBROUTINE  SKETCH  (NO , DCA , CA1 ) 

IMPLICIT  ;<eal*8  (a-h,o-z) 

COMMON  /BLKT8/G  (360,2) 

INTEGER  F  (101) .BLANK, DASH, TAB, DOT, STAR 

DATA  BLANK, DASH, TAB, DOT, STAR/'  1 | ‘ / 

1  FORMAT (2X, 'ave. 1 ,3X, 'S.D.  (deg  .)  ' , 1 01 A1 ) 

2  FORMAT  (lX,F6.i*,1X,F6.i»,1X,F5.1,101A1) 

3  FORMAT (//,T20, 'The  reference  unit  value  =  ' , E 1 2 . 6 , //) 
DO  100  1=1,101 

100  F  ( I ) =DASH 

WRITE  (6,1) F 
F  (1 )  =TAB 
DC  160  J=2 , 1 01 

160  F  (J) =BLANK 

UMAX=DABS  (G  (1  , 1)  )+DABS  (G  (1  , 1) ) 

REF=G  (1,1) 

DO  300  K=2 , NO 

TEM=DASS  (G  (K,1))+0ABS  (G(K,2)) 

I F (  UMAX  . LT .  TEM  )  UMAX=TEM 
I  F  (  REF  .LT.  DABS(G(K,1))  )  REF  =  DABS  (G  (K  ,  1 ) ) 

300  CONTINUE 

DO  500  L=1 , NO 

GL I NE=DABS  (G  (L , 1 ) /UMAX) 

L I NE=GL I NE*1 .D2+1 
IF  (  LINE  .LT.  2  )  LINE=2 
I F  (  LINE  .GT.  101  )  L I NE  =  1 01 
SL=DABS  (G  (L , 2) /UMAX) 

LS1=LINE+SL*1 .02 

IF  (  LS2  .LT.  2  )  LS 2=2 

IF  (  LSI  .GT.  101  )  LS 1 =1 01 

LS2=L I N  E — S  L  * 1 . D2 

IF  (  LS2  .LT.  2  )  LS2=2 

IF  (  LS2  .GT.  101  )  L S 2  =  1 01 

F  (LSI) =DOT 

F  (LS2) =DOT 

F  (L I NE) =STAR 
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CA=CA1+ (L  —  1 )  *0 C A 
GL I NE=G  (Lt1)/REF 
SL=G  (L , 2) /REF 
WRITE  (6,2) GL  I  NE  ,  SL , CA , F 
F  (L I NE) =BLANK 
F  (LSI ) =BLANK 
F  (LS2) “BLANK 
WRITE  (6,3) REF 
RETURN 
END 


C  *****1 

C 

C  *** 
C*** 

C*** 

(3  :'c  ;'c  ;'r 

c  *  *  * 
c  *** 

C*** 

C*  *>':**: 


CONVECTION  HEAT  TRANSFER  MODEL  OF  DIESEL  COMBUSTION 


SUBROUTINE  CONV  (OMEG , TG , TS) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKV2/R  (20)  ,QV  (20) 

COMMON  /BLK/RPM,RP,RL,RCS,CR,CM,HN,QC 
R  (1) =0.D0 
QV  (1 )  =0  .  DO 

CALL  PROPER (TG.RKG.VG) 

DO  100  1=2,11 

R  (I)  =.100*  (1-1)  *RP 
RE=2 .  D0*0MEG*R ( I )  *R ( I ) /VG 
QV  ( I ) = . 47D-1 *RKG/R ( I ) *RE** . 8D0* .700** .3300* (TG-TS) 
QV  ( I ) =QV  ( I ) /I .8805D2 
R (I )  -0. 1  DO*  (1-1) 

RETURN 


FLUID  PROPERTI ES 


SUBROUTINE  PROPER  (T , RKG , VG) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  RKGV(7)  ,VGV(7)  ,TV(7) 

DATA  RKGV/3.95DO, 5. 1300, 6.2500,7.32D0,8.37D0, 9-^600,1 .1101/ 

DATA  VGV/3. 790-1 ,6.61*0-1,9.90-1,1 .380+0,1 .82D+0,2.29D+0,3.1D+0/ 
DATA  TV/5 . D2.7-D2, 9-02,1 .103.1 -3D3.1-5D3.1 -803/ 

RKG=RKGV  (1 ) 

VG=VGV (1) 

IF  (  T  ,LT.  TV  (1 )  )  GOTO  900 
DO  100  1=2,7 

IF  (  T  ' .  LT .  TV  ( I )  )  GOTO  200 

CONTINUE 
RKG=RKGV (7) 

VG-VGV(7) 
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GOTO  900 

200  RKG=RKGV  (I  -1)  +  (RKGV  (I)  -RKGV  (I  -1) )  /  (TV  ( I )  -TV  ( I  -1)  )  *  (T-TV  (I  -1) ) 
VG=VGV  (I  -1)+  (VGV  (I)  -VGV  (I  -1)  )  /  (TV  ( I )  -TV  ( I  - 1 )  )  *  (T-TV  (I  -1)  ) 

900  RETURN 
END 


RADIATION  HEAT  TRANSFER  MODEL  OF  DIESEL  COM3UST I  ON 


C*******5 
C*** 

c*** 

C 

C  *** 

c*** 

C*** 

c*** 
c  * 

C*** 

SUBROUTINE  RANDO  (CAD) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR1/PI .PIST.PI4 

COMMON  /BLKR3/FVF.FVA,FVB,FVC,FAF , F A A , F AB , F AC , XHF , XHA , XHB , XHC , 
+  XCF , XCA , XCB , XCC 

COMMON  /BLKR4/RHODN , PH  I DN , ZDN , THE  I N , OMES . RNZ , VS , NFP 
COMMON  /BLKR5/RCB,ZCB,ZCC,ZCCB,RP 
COMMON  /BLKR7/REFL,REFTP,REFQ,REFW,REFWL,CVLWL 
COMMON  /BLKR8/OL,P,PE,T,TO,W 
COMMON  /BLKR9/AE  (20) 

COMMON  /BLKRA/RHOO , PH  I  0 , ZD , R I , STHE I N , CTHE I N 
COMMON  /BLKRD/TSCH.TSCW, TSP.TSB.ES 
.  DIMENSION  SAE  (10) 


C>»- 


c»> 

c»> 

c>» 

c»>- 
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I/O  FORMATS 


-<< 

<< 

<< 

<< 

-<< 


1 1 
12 
20 


21 


22 

23 


FORMAT  (6F  1  0 . 14 ,  / ,  F  1  0 . 4  ,  E  1  0 . 3 . 2 F  1  0 . 4 , 4  I  5 ,  / ,  2  I  5 . 3F  1  0 . 4  ,  /  ,  5F  1  0  .  L  , 
+  E10.3) 

FORMAT  (F10.4) 

FORMAT  (F10.4,E10.3,3F10.4,/,3 (10X.4F10.4,/)  ) 

FORMATC  Engine  speed  =  1  ,  F  5  •  0 ,  1  rpm'  ,/,  1  Combustion  bowl 
♦'radius  =  ',F6.4,'  cm',/,'  Compression  ratio  =  ',F4.1,/, 

+  '  Piston  radius  =  1 , F 6 . 2 , *  cm',/,'  Crank  sheft  length  =  ', 
+F6.2,'  cm',/,'  Connection  rod  length  =  '  ,  F  6  -  2 ,  ‘  cm',//, 

+'  Fuel  injection  rate  =  ',F6.4,'  cu.  cm/sec',/,'  Injection' 

+, '  nozzle  radius  =  ',F6.4,'  cm',/,'  Incline  angle  =  '.F4.1, 

+  '  deg',/,'  Swirl  ratio  =  ' ,F4.1,/,  '  Fuel  compostion 
+ '  C :  H :  Q  =  ',12,':', 12,':', 12,'  )',/,'  #  of  plumes  =  1 
FORMAT  (1 X ,' F 1 ame  temperature  -  AFT',/, 

+'  Surface  temperature  at',/,1  combustion' 

+, '  bowl  -  ' , F5 • 0 , ' K ' , / , 1  cylinder  head  -  ' , F5 .0, ' K ' ,/, 

+'  cylinder  wall  -  ' , F5.0, ' K ' ,/ , '  piston  surface  -  ', 
+F5.0, '«',//, '  Surface  emissivity  =  '  ,  F  4 . 2 ,  / ,  '  Length  and' 

+, '  wavelength  convertor  =  ' , F 1 0 . 0 , '  um/cm1,/) 

FORMATC  Engine  crank  angle  =  '  ,  F  6 . 1 ,  '  deg.  Piston', 

+  '  location  :  '  , F  6 . 4 , '  cm',/) 

FORMATC  Cylinder  pressure  =  '  ,  F  6 . 1  ,  '  KN/cm/cm '  ,  / ,  '  Distr 


( 


2,//) 
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+ 1 i but i on  paramters  of1,/,'  soot  volume  fraction  -  1  ,  E 1 0 . 2 , 
+3F10.2,/,'  burnt  fuel/air  ratio  - ' , F 1 0 . 4 , 3F 1 0 . 2 , / , 1  molar', 
+  '  fraction  of  H20  -  1  ,  F 1 0 . 4 , 3F 1  0 . 2  ,  / ,  20X ,  '  C02  -  '  ,  FI  0 .4, 3F1 0 .2) 

2U  FORMAT  (/,4X,  '  CA  =  '  ,  F7  -  2 ,  '  deg.,  P  =  ’,F6.2,'  atm, ' ,/,4X, 

+  'fv  =  ' , El  0 . 3. ' .  F/A  =  1 , F  8 . 5 , ' ,  XH20  =  ' , F8 . 5 , ‘ ,  XC02  =  ', 
+F8.5./) 

25  FORMAT  (4X, 'Spectral  AE  =  1 , 5E 1 0 . 3 ,/, 1 8X , 5E 1 0 . 3 ,/, 4X , 1  Apparent ' 
+  , '  Emissivity  =  ■ 1 , F 8 . 5 » /) 

26  FORMAT  (IX, 'Detector1  's  location  -  ( 1 , F A  .  1 , ' , 1 , F A . 0 ,  '  ,  1  ,  F 4  .  1  , 


+  ')  ') 

C»> - <« 

C»>  «< 

C»>  INPUT  :  engine  speed  (RPM)  ,  radius  of  combustion  chamber  <<< 

C»>  (RCB)  ,  compression  ratio  (CR)  ,  piston  radius  (RP)  <<< 

C»>  length  of  engine  crank  sheft  and  connection  rod  <<< 

C»>  (RCS.RL),  fuel  injection  rate  (FIR),  injection  <<< 

C»>  nozzle  radius  (RNZ)  ,  fuel  injection  angle  (THE1N)  «< 

C»>  ,  swirl  ratio  (OMES)  ,  fuel  composition  (CM.HN.OI)  «< 

C»>  ,  §  of  fuel  plumes  (NFP)  ,  detector's  location  «< 

C»>  (RHODN ,  PH  I  DN  ,  ZDN)  ,  surface  temperatures  (TSB,  <<< 

C>»  TSCH , TSCW.TSP)  ,  surface  emissivity  (ES)  and  <« 

C»>  unit  convertion  factor  (CVLWL)  .  «< 

C»>  «< 

C»> - <« 


READ (5,10)  RPM, RCB, CR.RP, RCS.RL 
+  ,F  IR, RNZ, THE  I N , OMES , MC , NH , 10, NFP 

+  ,NTHE, NOP, RHODN, PHIDN, ZDN 

+  , TSB, TSCH, TSCW.TSP.ES, CVLWL 

WRITE (6,20)  RPM, RCB, CR.RP, RCS.RL 
+  ,FIR, RNZ, THEIN, OMES, MC.NH, 10, NFP 

WRITE (6, 21)  TSB, TSCH, TSCW.TSP.ES, CVLWL 
CM=MC 
HN=NH 
01*10 
REFL=RP 

C»> 

C»>  Normalize  variables  by  reference  unit. 

C»> 

RP=RP/REF  L 

RCB=RCB/REFL 

RCS=RCS/REFL 

RL=RL/REFL 

RNZ=RNZ/REFL 

THE IN-THE IN*PI/.l8D3 

RHODN=RHODN/REFL 

PH  I DN=PH I DN*P I / . 1 8D3 

ZDN=ZDN/REFL 

ZCB=2.00*RP*RP*RCS/RCB/RCB/  (CR- 1 .DO) 

F I R*F I R/REFL/REFL/REFL 

VS*F IR/PI/RNZ/RNZ/NFP 

0MES=0MES*RPM*PI/3.D1 

TS8=TSB/REFTP 

TSCH*TSCH/REFTP 

TSCW-TSCW/REFTP 

TSP*TSP/REFTP 
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STHE I N=DS I N (THE  I N) 
CTHE ! N=DC0S (THE  IN) 

C»> - 

C»> 


CHECK  if  data  ends  -  READ  engine  crank  angle  (deg) 
IF  CA  >  0,  go  on  computing 
IF  CA  <=  0,  stop. 


READ  (5.11) CAD 

IF  (  CAD  .LE.  0 . DO  )  GOTO  900 
CA=CAD*P I / . 1 8D  3 

CW=DSQRT  (1 . DO-RCS*RCS/RL/RL*DS I N (CA) **2) 
ZCC=2.D0*RCS- (RL*CW-RCS*DCOS (CA)+RCS-RL) 
ZCCB=ZCC+ZCB 


READ  cylinder  pressure  (P)  ,  distribution  parameters 

of  soot  volume  fraction  (FVF , FVA , F VB , F VC) ,  burnt 
fuel/air  ratio  (F AF , F AA, F AB , F AC)  ,  water  vapor  ( 
(XHF , XH A , XHB , XHC) ,  and  carbon  dioxide  (XCF.XCA, 
XCB.XCC) 


READ  (5.12) P,FVF,FVA,FVB,FVC,FAF,FAA,FAB,FAC 
,  XHF,  XHA, XHB, XHC, XCF.XCA, XCB.XCC 
WRITE  (6,23) P,FVF,FVA,FVB,FVC,FAF,FAA,FAB,FAC 
, XHF, XHA, XHB, XHC, XCF.XCA, XCB.XCC 
REFP=1 .0332300 
P=P/REFP 


COMPUTE  time-resolved  radiation  heat  flux  from 
in-cylinder  species  to  detector  (RHF) 

WRITE  engine  crank  angle  (CA) ,  mean  soot  volume  fraction 
(FVM)  ,  mean  burnt  fuel/air  ratio  (FAM)  ,  mean  molar 
fraction  of  water  vapor  (XHM) ,  mean  molar  fraction 
of  carbon  dioxide  (XCM) ,  and  apparent  emissivity 
E 


DO  800  1=1,11 

RH0DN=.1D0*  (1-1) 

ZDN=0 . DO 
DO  800  J=1 ,1 

PH  I DN= .501  *  (J-1 ) 

CALL  XMEAN (FVM, FAM, XHM, XCM) 

IF  (  I  .LE.  1  )  WRITE (6, 2L)  CAD, P, FVM, FAM, XHM, XCM 
WRITE  (6,26)RHODN,PHIDN,ZDN 
PH  I DN=PH I DN*P I / . 1 8D  3 
AE  ( I ) =RHF  (NTHE , NOP , CM, HN , 0 1  ,SAE) 

800  WRITE  (6,25) SAE.AE (I) 

900  RETURN 
END 

C  *  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  i 
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C**ft 

C*** 

C  *  *  * s': s'; s': s':  s';  *  s':  s';  V;  V;  s';  s';  s'; 
c»> - 

c»> 

c»> 

c»> 


COMMON  FUNCTIONS 


C3C>r2C5C5C;C5C5CT 


*  J#  »'«  A  A  A  A  A  A  A  A  A  A  A  ** 

»  It  A  It  It  It  It  It  It  It  It  It  It  It  II 


PERMUTATION 


1 
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DOUBLE  PRECISION  FUNCTION  PER  (M, N) 

IMPLICIT  REAL*8  (A-H.C-2) 

PER=1 .DO 

IF  (  M  .LT.  N  )  GOTO  200 

IF  (  (M  .EQ.  0)  .OR.  (N  .EQ.  0)  )  GOTO  900 
GOTO  300 
WRITE  (6,201) 

FORMAT (IX, 1  ERROR  in  PERMUTATION  function') 
STOP 

Ml =M-N+1 

DO  310  I =M1 , M 

PER=PER*I 

I F  (  PEP.  .GE.  1  .D70)  GOTO  200 

CONTINUE 

RETURN 

END 


RADIATION  INTEGRAL  FUNCTION 


1983  “  Chang  and  Rhee 


DOUBLE  PRECISION  FUNCTION  RlF(M,A,X) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  Y  (10) 

FORMAT  (IX, 'ERROR  in  RIF  f unc t i on . ' , // , 1 X , 1 1 0 , 2 E 1 5 . 4) 

R I F=0 . DO 

N=0 

IF  (  (  M  .LT.  0  )  .OR.  (  A  .LT.  0.D0  )  )  GOTO  200 

IF  (  M  .GT.  9  )  GOTO  200 

IF  (  X  )  200,400,300 

WRITE (6, 1) M.A.X 

STOP 

Ml *M+1 

DO  310  I  1=1 , Ml 
1=11-1 

Y (1+1) =PER  (M,  I) 

N=N+1 
SUM=0.D0 
DO  3^0  1 1=1 , Ml 
1*11-1 

SUM=SUM+Y  ( I +1 ) *X** (M- I )  / (A+N)  ** ( I +1 ) 

SUM=SUM*DEXP  (-  (A+N) *X) 

R I F*R I F+SUM 

IF  (  DABS  (SUM) -.ID-10  )  900,900,320 
N-.5D6**  (1 .DO/  (M+1 ) ) -A 


*** 
*  *  s'; 

3,?:,C3'C3,C3,C:,C3%3,CJ,C 

- <<< 

«< 

«< 

«< 

- <<< 


«< 

«< 

«< 

«< 

«< 

«< 

<« 
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■  %  P  s.  ' 

*  rnUktm 


,J  -jf  S 
le-VCtc. 


.  -  .  ••  .v  ' 


SUM-0  .DO 
DO  420  1=1 ,N 

420  SUM-SUM+1 .DO/  (A+l)  ** (M+1) 

R I F  =PER (M.M) * (SUM+.5D0* (1 .00/ (A+N+1)**M+1 .DO/ (A+N)  **M) /M) 
900  RETURN 
END 


C»> - «< 

C>»  <« 

C»>  ERROR  (PROBABILITY  INTEGRAL)  FUNCTION  «< 

C»>  «< 

C»> - «< 


DOUBLE  PRECISION  FUNCTION  ERRF  (X) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR1 /P I , P I  ST , P I  4 
COMMON  /BLKR2/CFRE  (40) .ERFC.ERFCL 
Y=DABS  (X) 

I F  (  Y  .GT.  1 .0-15  )  GOTO  200 

ERRF=X 

ERF  C=1 .DO 

ERF  CL— X-X*X/2 .  DO 

IF  (  Y  .LE.  .990-15  )  GOTO  900 

ERRF 1=ERRF 

200  IF  (  Y  .GE.  2.3500  )  GOTO  400 
TERM=2 . 00*X/P I  ST 
ERRF=TERM 
A=0 .  DO 

300  A-A+1 . DO 

.  TERM— TERM*X*X*(2.D0*A-1  .DO)  /A/  (2.D0*A+1  .DO) 
ERRF-ERRF+TERM 

IF  (DABS  (TERM)  .GT.  .ID-10  )  GOTO  300 
IF  (  Y  .GE.  1.D-15  )  GOTO  320 

ERRF-(ERRF*(1.D0-Y*1.D15)+ERRF1*(Y*1  .D15-.9900))*1  .02 
320  ERFC-X/Y-ERRF 

ERFCL=OLOG (DABS (ERFC) ) 

IF  (  Y  .LE.  2.34D0  )  GOTO  900 
ERRF1=ERRF 
400  B=X*X-.5D0 

IF  (  Y  .GT.  8. DO  )  GOTO  800 
TERM=DEXP  (-X*X) /X/PIST 
ERF  C=TERM 
I  B=B 

DO  600  IA=2, IB 

IF  (  DABS  (TERM)  .LT.  1.D-7)  GOTO  700 
TERM— TERM*  (2*1  A-3)  *-5D0/X/X 
600  ERFC-ERFC+TERM 
700  ERRF=X/Y-ERFC 

IF  (  Y  .GE.  2.3500  )  GOTO  720 
ERRF- (ERRF* (2. 35DO-Y)+ERRF1*(Y-2. 3400) ) *1 .02 
720  ERFCL=DLOG  (OABS  (ERFC) ) 

IF  (  Y  .LE.  7.9900  )  GOTO  900 
ERFC1-ERFC 
800  ERRF-Y/X 

ERF  CL— X*X-DLOG  (P I  ST) 

I B=X*X- .500 
TERM-1 .DO/X 
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SUM=TERM 
DO  820  1=2, IB 

IF  (  DABS  (TERM)  .LT.  .ID-10  )  GOTO  840 
TERM=TERM* (2*1-3) *.5D0/X/X 
820  SUM=SUM+TERM 
840  ERFCL  =  ERFCL+DLOG (DABS  (SUM) ) 

£RFC=1 .0-75 

IF  (  DABS  (ERFCL)  .LT.  .1  72703  )  ERF  C=DEXP  (ERFCL) 
ERFC=Y/X*ERFC 

IF  (  Y  .GE.  8.00  )  GOTO  900 
ERF  C= (ERF  C* (8 . DO-Y) +ERFC1 * (Y-7 . 9900) ) *1  .02 
900  RETURN 
END 


C»> - <<< 

C»>  <<< 

C»>  COEFFICIENTS  OF  FRE  EXPANSION  «< 

C»>  <<< 

c»> - «< 


SUBROUTINE  I CFRE  (NTERMS) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR1/PI .PIST.PI4 
COMMON  /BLKR2/CFRE  (40) ,ERFC, ERFCL 
DIMENSION  G  (2,41) 

G  (1  ,1)=PIST*.5D0 
00  300  N=1, NTERMS 


A=2 . D0*N~3 . DO 
N 1 =N+1 

G  (1  , N1 )  =0 . DO 
DO  100  I =1 , N1 
G  (2 ,  I )  =G  (1 ,  I ) 

100  CONTINUE 

110  IF  (N-1)  900,210,120 

120  G (1 ,N) =4 . DO*A*G  (2, N-1) /  (A+2.D0) *PI/4.D0 

G  (1 , 1)  =  (G  (2 , 2)  +A*G  (2 , 1 ) )  /  (N-1 )  /  (A+2  .DO)  *P  I  /4 . DO 
130  IF  (N-2)  210,210,140 

140  N2=N-1 

DO  200  1=2, N2 

G  (1  ,  I)  =G  (2,  1+1)  *1*  (2*1-1) 

G  (1 ,  I)  =  (G  (1  ,  l)  +  (2.00*  (2 .  D0*A+1  .DO)  *  (I  -1)  +A)  *G  (2,  I) )  /N2 
200  G  (1 ,  l)  =  (G  (1  ,  I)  /  (A+2 . DO)  +4 . DO*A*G  (2,1-1)/  (A+2. DO) )  *PI/4.D0 

210  CFRE  (N)  =G  (1  , 1 ) 

300  CONTINUE 
900  RETURN 
END 

C»> - <<< 

c»>  <<< 

c»>  INVERSE  ERROR  FUNCTION  «< 

C»>  <« 

c»> - <<< 

DOUBLE  PRECISION  FUNCTION  FRE  (Y , YC , YL) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR1/PI .PIST.PI4 
COMMON  /BLKR2/CFRE  (40)  ,ERFC, ERFCL 
FRE=0 . DO 


X=DABS  (Y) 

XC=OABS  (YC) 

I F  (  X  .GT.  1 .DO  )  GOTO  150 
IF  (X  .LE.  .900  )  GOTO  200 
IF  (X  .LE.  . 99953DO)  GOTO  300 
IF  (  XC  .GE.  1 .0-75  )  GOTO  400 
FRE=. 26282802 
GOTO  900 

150  WRITE(6,15D 

151  FORMAT (IX. 1  ERROR  in  FRE  function.1) 

STOP 

200  A=X 

FRE=FRE+A*CFRE  (1) 

210  DO  230  1=2,40 
A=A*X*X 

FRE=FRE+A*CFRE  (I) 

IF  (  A*CFRE(I)  . LT.  1.0-7  )  GOTO  900 
230  CONTINUE 
240  GOTO  900 
300  DO  320  l=i ,100 

A=PIST*.5D0*DEXP (FRE*FRE) * (X-ERRF (FRE) ) 
FRE=FRE+A 

IF  (  DABS  (A)  .LT.  1.0-7  )  GOTO  9°0 
320  CONTINUE 
GOTO  900 
400  FRE=2 . DO 

DLP I ST=DLOG  (P I  ST) 

’  DLXC=DLOG  (XC) 

IF  (  XC  .LE.  1.D-75  )  0LXC=YL 
DO  500  J=1 , 20 

IF  (  A  .LT.  FRE*rRE- .5DO  )  A=FRE*FRE- . 500 
TERM=1 .DO/FRE 
SUM=TERM 
I  A=A 

DO  430  1=2, I  A 

TERM-TERM*  (2*1-3)  * . 5D0/FRE/FRE 
SUM=SUM+TERM 

IF  (  DABS  (TERM/SUM)  .LT.  .ID-10  )  GOTO  440 
430  CONTINUE 

440  FRE1=FRE 

TERM=DLOG (SUM) -DLXC-DLPI ST 
IF  (  TERM  .LT.  O.DO  )  GOTO  900 
FRE=DSQRT  (TERM) 

IF  (DABS  (FRE-FRE1)  .LE.  .10-10  )  GOTO  900 
500  CONTINUE 
900  FRE=Y/X*FRE 
RETURN 
END 


C»> - <<< 

C»>  <« 

C»>  GAUSS  NUMERICAL  INTEGRATION  «< 

C»>  «< 

c»>  «< 

c»>  <<< 

c»> - <<< 
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DOUBLE  PRECISION  FUNCTION  GNI(XU,XL,Y) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  C  (5)  ,Y(5) 

DATA  C/.1184634425DO, . 23931 43352DO , . 2844444444DO , . 23931 43352D0 , 
i-  .118463442500/ 

GN I =0 . DO 
DO  100  1=1,5 
GNI=GNI+C  (I)  *Y  (I) 

GN I =GN I  *  (XU-XL) 

RETURN 
END 


X  positions  of  5  points  GN I 


SUBROUTINE  GN I  X  (XU , XL , X) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  U  (5)  ,X  (5) 

DATA  U/-. 906179846000, -.5384693102DO.O. DO, .53846931 02D0 , 
+  .906179846000/ 

DO  100  1=1 ,5 

X  ( I ) = . 5D0*  (XU+XL+  (XU-XL)  *U ( I ) ) 

RETURN 

END 


Adiabatic  Flame  Temperature  coefficients 


SUBROUTINE  AFTC  (P , CM, HN , 0  I ) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR7/REFL,REFTP,REFQ,REFW,REFWL,CVLWL 
COMMON  /BLKRC/TA,TB,TC 

FORMAT  ( 1  ERROR  in  AFTC  -  P  =  '.E12.4,'  atm.  Fuel  =  \3E12.4) 
IF  (  (CM  .GE.  0  .DO)  .AND.  (HN  .GE.  0.D0)  .AND.  (01  .GE.  0.D0) 

+  .AND.  (P  .GE.  0 .DO)  )  GOTO  100 
WRITE  (6,1) P , CM, HN , 0  I 
STOP 

IF  (  01  .LE.  0 . DO  )  GOTO  200 
A 1 cohol  f ami  1 y 

IF  (  (CM  .LT.  .1901)  .OR.  (CM  .GT.  . 8 1 D 1 )  )  GOTO  90 
A1=.21535D4+.602D2*DLOG (CM-. 1 51 01) 

A2=-.5D0+. 5201*0 LOG  (CM+.1100) 

B1  =  .641  DO- . 3D-1*DLOG (CM+.34D1) 

B2=.4D-3+.66D-2*DLOG  (CM+.2D-1) 

C1=.193DO-.41D-1*DLOG (CM+.35D1) 

C2—  .2D-3+.66D-2*DLOG  (CM+.9D-1) 

GOTO  900 

Hydrogen 
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200  IF  (  CM  .GT.  0 . DO  )  GOTO  300 
A 1  = .  2  4 1  2  2  D  4 
A2=. 1 79D2 
B1  =  .  1*7300 
B2= .210-1 


Cl = . 21 7D-1 
C2=. 1 91 D- 1 
GOTO  900 

C>>> 

C»>  Acetylene 

C»> 

300  IF  (  (CM  .NE.  2. DO)  .AND.  (HN  .NE.  2. DO)  )  GOTO  400 
A1=. 2562904 
A2= . 423D2 
B1=.39D0 
B2=.426D-1 
Cl =- . 81 9D-1 
C2= . 379D- 1 
GOTO  900 

C»> 

C»>  Paraffin  family 

C»> 

400  IF  (  HN  . LT .  2.D0*CM+1.00  )  GOTO  500 
IF  (.CM  .GT.  .  1 61 D2  )  GOTO  90 
A1  = . 2295D4+ . 950 1 *DLOG (CM-. 9900) 

A2-.13202+.7D0*DLOG (CM-.97DO) 

B1 ». 54400-. 19D-2*DLOG  (CM-. 9900) 

B2=.18D-1  +  .9D-3*DI_0G  (CM-.96D0) 

C1-.7D-1-.3D-2*DLOG (CM-.98DO) 

C2=.1 690-1+. 8D-3*OLOG (CM-. 9700) 

GOTO  900 

500  IF  (  HN  .LT.  2.D0*CM-1 .DO)  GOTO  600 

C»> 

C»>  Olefin  family 

C»> 

IF  (  (CM  .LT.  .1901)  .OR.  (CM  .GT.  .161D2)  )  GOTO  90 
A1*.2377D4-.179D2*DLOG (CM-.I85DI) 
A2=.19902-.14D1*DLOG (CM-.189DI) 

B1=.512DO-.85D-2*DLOG (CM-.I89DI) 
B2=.250-1-.14D-2*DLOG (CM-.186D1) 
C1=.350-1+.83D-2*DLOG (CM-.189D1) 
C2-.23D-1-.130-2*DLOG (CM-.18601) 

GOTO  900 

C»> 

C»>  Aromatic  family 

C»> 

600  IF  (  (CM  .LT.  . 59D1 )  .OR.  (CM  .GT.  .22102)  )  GOTO  90 
A1-. 237804-. 153D2*DLOG (CM-.536DI) 

A2*.21D2- .16D1*DL0G (CM-. 52901) 

B1«.5l8DO+.6D-2*DLOG (CM- .55101) 

B2-.  270-1 -.2D-2*DL0G (CM-. 52401) 

Cl =. 340-1+. 8D-2*DLOG (CM-. 53801) 

C2* . 250-1  - . 1 6D- 2*0 LOG (CM-. 52701) 


o 


Estimated  AFT  Formulas 
X=OLOG  (P) 

TA=  (A1+A2*X)  /REFTP 

TB=B1+B2*X 

TC=C1+C2*X 

RETURN 

END 


Adiabatic  Flame  Temperature  estimate 


DOUBLE  PRECISION  FUNCTION  AFT  (FA) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKRC/T A ,TB , TC 

IF  (  FA  .GE.  0 . DO  )  GOTO  100 

WRITE  (6,1) FA 

FORMAT  ( 1  Error  in  AFT  :  F.'A=' ,  El  2 . 3) 

STOP 

IF  (  (FA  .LE.  .400)  .OR.  (FA  .GE.  1.6DO)  )  GOTO  200 
F  ALOG=DLOG (FA) 

IF  (  FA  .GT.  1  . DO)  FALOG=DLOG  (2. DO-FA) 

AFT=TA*  (1 .DO+TB*FALOG+TC*FALOG*FALOG) 

RETURN 

AFT=.25D0 

RETURN 

END 


Mean  distribution  properties 


SUBROUTINE  XMEAN  (FVM, F AM, XHM, XCM) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR1/PI , P I  ST , P I  4 

COMMON  /BLKR3/FVF,FVA,FVB,FVC,FAF , F AA , F AB , F AC , XHF , XHA , XHB , XHC , 
i-  XCF  ,  XCA ,  XCB ,  XCC 

COMMON  /BLKR4/RHODN , PH  I  ON , ZDN .THE  I N , OMES , RNZ , VS , NF P 
COMMON  /BLKR5/RCB,ZCB,ZCC,ZCCB,RP 
COMMON  /BLKRA/RHOD , PH  I D , ZD , R I , STHE I N , CTHE I N 
DIMENSION  RX(5)  ,WX  (5)  ,ZX(5) 

DIMENSION  F  V I  (5)  ,  F VJ  (5)  ,  FVK  (5)  ,  F  A  I  (5)  ,FAJ(5)  ,FAK(5) 

DIMENSION  XHI  (5)  ,XHJ  (5)  ,XHK  (5)  ,  XC  I  (5)  ,XCJ  (5)  ,XCK  (5) 
VOL»PI*(ZCC*RP*RP+ZCB*RCB*RCB) /NFP 
FVM=0 . DO 
F  AM=0 . DO 
XHM-O.DO 
XCM-0.00 
DO  20  L-1,5 
FVI  (L) =0 . DO 
FVJ  (L) =0 . DO 
FVK  (L)  "0 .  DO 
FA  I  (L)  *0 .  DO 
FAJ  (L) *0 . DO 


**  •"“*  t*  .*•"  •  V,  -*■  ' 
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FAK  (L)  =0 . DO 
XHI  (L) =0 . DO 
XHJ  (L) =0 .DO 
XHK  (L) =0 .DO 
XCI  (L) =0 .00 
XCJ  (L) =0 . DO 
20  XCK  (L)  =0 .DO 

RU=RP 
WU=P I / N  F  P 
ZU=0.00 
ZL=-ZCC 

1*0  CALL  GNIX  (ZU.ZL, ZX) 

DO  300  1=1,5 
Z= ZX  ( I ) 

CALL  GNIX(WU.O.DO.WX) 

DO  200  J=1 ,5 
W=WX  (J) 

SW=DS I N (W) 

CW=DCOS (W) 

CALL  GNIX  (RU.O.DO.RX) 

DO  100  K-1,5 

R=DSQRT  (RX  (K) ) 

X=R*CTHE I N*CW-Z*STHE I N 
Y=R*SW 

RN=DSQRT  (X*X+Y*Y) 

WN=DATAN2  (Y ,  X) 

ZN=R*STHE IN*CW+Z*CTHE IN 

G=-  ( F  V  A  *  R  N+ F  VB  *  WN  *  WN+ F  V  C  *  Z  N  *  Z  N ) 

FVK  (K)  =FVF*DEXP  (G) 

G=-  (FAA*RN+FAB*WN*WN+FAC*ZN*ZN) 

FAK(K)=FAF*DEXP(G) 

G=-  (XHA*RN+XHB*WN*WN+XHC*ZN*ZN) 

XHK  (K) =XHF*DEXP  (G) 

G=-  (XCA*RN+XCB*WN*WN+XCC*ZN*ZN) 

100  XCK  (K)=XCF*DEXP  (G) 

F  VJ ( J) =GN I  (ZU.ZL.FVK) 

FAJ  (J) =GN I  (ZU.ZL, FAK) 

XHJ  (J) =GN I  (ZU.ZL, XHK) 

200  XCJ  (J) =GN I  (ZU.ZL, XCK) 

FVI  (I) =GNI  (WU.WL.FVJ) 

F  A  I  ( I ) =GN I  (WU , WL , F  AJ) 

XHI  ( I ) =GN I  (WU.WL.XHJ) 

300  XCI  (l)*GNI  (WU,WL,XCJ) 

FVM=FVM+GN I  (RU , 0 . DO , F V I ) /VOL 
FAM-FAM+GNI  (RU.O.DO.FAI) /VOL 
XHM=XHM+GN I  (RU , 0 . DO , XH I ) /VOL 
400  XCM=XCM+GN I  (RU , 0 . DO , XC I ) /VOL 

IF  (  DABS  (ZL+ZCCB)  .LT.  .10-10  )  GOTO  900 
ZU=-ZCC 
ZL— ZCCB 
RU=RCB*RCB 
GOTO  40 
900  RETURN 
END 

C  A  ft  *  ft  ft  ft  A  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  ft  i 


/•'j  '•'.’-".V.".-' "-“.V. „V •.•fly’ /». ’’ .V." 


RADIATION  HEAT  FLUX 


AAA 


DOUBLE  PRECISION  FUNCTION  RHF  (NTHE , NOP , CM, HN , 0  I ,SAE) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR1/P! ,PIST,Pl4 

COMMON  /BLKR4/RHODN,PHIDN,ZDN,THEIN,OMES,RNZ,VS,NFP 
COMMON  /BLKR8/OL,P,PE,T,TO,W 
COMMON  /BLKRA/RHOD,PHID,ZD,RI , STHE I N , CTHE I N 
COMMON  /BLKRC/TA.TB.TC 

COMMON  /BLKRE/THEN , ZETN , STHEN , CTHEN , SZETN , CZETN 
COMMON  /BLKRF/QRS  (10)  ,RHF RNR , QRSNR  (10)  ,RHFRP,QRSP  (10) 

DIMENSION  CTHEX  (5)  ,ZETX  (5)  ,SAE  (10) 

DIMENSION  QTZ(5)  ,QSTZ(5,10)  ,  QT  (5)  .  GQS  (1  0)  ,QST(5.10) 

DIMENSION  QTZNR  (5)  ,  QSTZNR  (5 , 1  0)  ,  QTNR  (5)  ,GQSNR(10)  ,QSTNR(5,10) 
DIMENSION  QTZP(5)  ,QSTZP(5,10)  ,QTP(5)  .GQSP(IO)  ,QSTP(5,10) 
DIMENSION  SRHF  (10)  ,  SRHFNR  (10)  ,  SRHFP  (10) 

1  FORMAT  (IX,  1  QN  =  1  , 6E1 0 . 3 . / .  1 8X , 5E 1  0 . 3 , /,  1 X ,  1  QSW  =  '.6E10.3, 
+/.18X.5E10.3./.1X, ‘QSP  =  ' ,6Ei0.3,/,i8x,5Eio.3) 

RHF=0 . DO 
RHFNR=0 . DO 
RHF  P=0 . DO 
DO  100  1=1,10 
SRHF  (I) =0.D0 
SRHFNR  (I) =0. DO 
100  SRHFP  ( I ) =0 . DO 
C»> 

C>»  Compute  flame  temperature  parameters. 

C»> 

CALL  AFTC (P , CM, HN , 0  I ) 

IF  (  DABS  (PH  I DN)  .LE.  PI/NFP  )  GOTO  120 
JFP=DABS  (PH  I DN) *NFP/2 . DO/P I  +  .5D0 
PH  I FP=P l*2.D0*JFP/NFP 
IF  (  PH  I DN  .LT.  0 . DO  )  PH l F  P=-PH I FP 
PH  I DN=PH I DN+PH I FP 

C»> 

C»>  Compute  heat  flux  in  specified  solid  angle. 

C»> 

120  CTHEL-+1.D0 

DCTHE=-2 . DO/NTH E 

C»> 

C»>  Integrate  heat  flux  in  sub- i  n  terva  1  s  of  THE  domain. 

C»> 

DO  600  J=1 , NTHE 

CTHEU=CTHEL+DCTHE 

CALL  GNIX  (CTHEU.CTHEL, CTHEX) 

DO  500  K-1,5 
QT  (K) =0.00 
QTNR  (K) =0  .DO 
QTP  (K) =0.00 
DO  140  K1 =1 , 10 
QST  (K , K1 ) =0 . DO 
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140 


QSTNR  (K,K1) =0.D0 
QSTP  (K , K1 ) =0 . DO 
2ETL=-P I 
CTHEN-CTHEX  (K) 

STHEN=DSQRT (1 . D0-CTHEN*CTHEN) 

IF  (  STHEN  .LE.  O.DO  )  GOTO  500 
THEN=DATAN2  (STHEN , CTHEN) 

NZET=1+STHEN/.25D0 

IF  (  (RHODN  .LT.  1.D0)  .AND.  (STHEN  .GT.  -75DO)  )  NZET=2 
IF  (  RHODN  .LT.  .ID-10  )  NZET-NZET+NZET 
DZET=2 .  DO*P I /NZET 

C»> 

C»>  Integrate  heat  flux  in  sub- i  nterva  1  s  of  ZET  domain. 

C»> 

DO  400  L=1 , NZET 
ZETU=ZETL+DZET 

IF  (  DABS  (ZETU)  .LT.  .ID-10  )  ZETU=0.D0 
CALL  GN I  X  (ZETU , ZETL , ZETX) 

DO  300  M-1,5 
ZETN=ZETX (M) 

IF  (  DABS  (ZETN)  .LE.  .ID-10  )  ZETN-O.DO 
SZETN-DS I N (ZETN) 

CZETN-DCQS (ZETN) 

0THEN-THEN*.18D3/PI 
DZETN-ZETN*.l8D3/PI 
PDN=PH I DN 

QTZ  (M) =RHFR  (NOP , PDN , REOP) 

QTZNR(M) =RHFRNR 
*  QTZP (M) -RHFRP 

DO  200  MI-1,10 

QSTZ  (M.MI)-QRS(MI) 

QSTZNR  (M,M1) -QRSNR (Ml) 

200  QSTZP (M,M1)  =QRSP (Ml) 

300  CONTINUE 

IF  (  REOP  .LT.  .10-10  )  GOTO  400 
QT  (K) -QT  (K) +GN I  (ZETU, ZETL , QTZ) /P I 
QTNR  (K)  -QTNR  (K)+GNI  (ZETU ,  ZETL  ,  QTZNR)  /PI 
QTP  (K) -QTP  (K) +GN I  (ZETU, ZETL , QTZP) /P I 
DO  320  K 1 = 1 ,10 
DO  310  K2=1 ,5 

GQS (K2) =QSTZ  (K2,K1) 

GQSNR (K2)  =QSTZNR  (K2,K1) 

310  GQSP  (K2) -QSTZP  (K2.K1) 

QST(K.KI)  =QST(K,K1)+GNI  (ZETU ,  ZETL  ,  GQS)  /PI 
QSTNR  (K,K1) -QSTNR  (K,K1)+GNI  (ZETU , ZETL , GQSNR) /PI 
320  QSTP (K,K1) -QSTP (K,K1)+GNI  (ZETU , ZETL , GQSP) /P I 

400  ZETL-ZETU 

500  CONTINUE 

RHF-RHF-GN I (CTHEU , CTHEL , QT) 

RHFNR-RHFNR-GN I (CTHEU , CTHEL , QTNR) 

RHFP-RHFP-GN I  (CTHEU , CTHE L , QTP) 

DO  520  LI-1 ,10 
DO  510  L2-1.5 

GQS  (L2) -QST  (L2 , LI) 

GQSNR  (L2) -QSTNR  (L2, LI) 


510 


GQSP  (L2) =QSTP  (L2, LI) 

SRHF  (LI) =SRHF  (LI) -GNI  (CTHEU , CTHEL , GQS) 

SRHFNR  (L1)=SRHFNR  (LI) -GNI  (CTHEU, CTHEL , GQSNR) 
520  SRHFP  (L1)=SRHFP  (LI) -GNI  (CTHEU , CTHEL , GQSP) 

600  CTHEL=CTHEU 

WRITE (6,1) RHF , SRHF, RHFNR, SRHFNR, RHFP, SRHFP 
Rl F 1=1 .DO 
DO  700  1=1,10 
X=.7D0*I 

R I F2=R IF (3 , 0 . DO , X)  * . 1 5D2/P  I  U 
IF  (  I  .GE.  10  )  R I F 2=0 . DO 
SAE  ( I )  =SRHF  (I)  /  (Rl  F1-RI  F2) 

700  R I F 1 =R I F  2 


RETURN 

END 

C»> - «< 

C»>  <<< 

C»>  RANGE  of  an  OPTICAL  PATH  «< 

C»>  «< 

c»> - <« 


SUBROUTINE  ROP  (REOP , TDS , TES , ODS , OE S , N) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR1 /P I , P I  ST , P I k 

COMMON  /BLKRi*/RHODN,PHIDN,ZDN,THE  IN.OMES.RNZ,  VS.NFP 

COMMON  /BLKR5/RCB,ZCB,ZCC,ZCCB,RP 

COMMON  /BLKRA/RHOD , PH  I D , ZD , R I , STHE I N , CTHE I N 

COMMON  /BLKRD/TSCH,TSCW,TSP,TSB,ES 

COMMON  /BLKRE/THEN , ZETN , STHEN , CTHEN , SZETN , CZETN 

RE0P=0 . DO 

TDS=0 . DO 

TES=0 .DO 

0DS=0 . DO 

0ES=0 .DO 

PH  I N=0  .  DO 

A=DCOS  (PH  I DN+ZETN) 

C»> 

C»>  Check  if  detector  locates  on  combustion  surface  and 
C>>>  optical  ray  shoots  into  combustion  chamber. 

C»> 

IF  (  (-ZDN-ZCC  .GE.  -.ID-10)  .AND.  (DABS (ZCC)  .GT.  .ID-10)  ) 

+  GOTO  200 

IF  (  (DABS  (ZCC)  .LT.  .ID-10)  .AND.  (  RHOON  .GT.  RCB)  )  GOTO  200 

IF  (  (DABS  (ZCC)  .LT.  .ID-10)  .AND.  (RHODN  .LE.  RCB)  )  GOTO  210 

IF  (  (ZDN  .GT.  .10-10)  .OR.  (RHODN-RP  .GT.  .ID-10)  )  GOTO  900 
IF  (  (DABS  (ZDN)  .LT.  .ID-10)  .AND.  (THEN  .GE.  .5D0*P|)  )  GOTO  900 
IF  (  (A  .GE.  0 . DO)  .AND.  (DABS  (RHODN-RP)  .LT.  .ID-10)  )  GOTO  900 

IF  (  (DABS  (ZDN+ZCC)  . LT .  .10-10)  .AND.  (RCB-RHODN  .GT.  SAMLL)  ) 

+  GOTO  900 

IF  (  (DABS  (ZDN+ZCC)  .LT.  .10-10)  .AND.  (THEN  .LT.  .5D0*PI)  ) 

+  GOTO  900 
GOTO  210 

200  IF  (  (-ZDN-ZCCB  .GT.  .ID-10)  .OR.  (RHODN-RCB  .GT.  .ID-10)  ) 

+  GOTO  900 

IF  (  (DABS  (ZDN+ZCCB)  .LT.  .ID-10)  .AND.  (THEN  .LE.  .5D0*PI)  ) 

+  GOTO  900 


c»> 

c»> 

c>» 

210 


220 

C»> 

c»> 

c»> 

300 


310 

320 


330 

340 

400 


410 


420 


430 

440 

450 


460 

510 


IF  (  {DABS  (RHODN-RCB)  .LT.  .10-10)  .AND.  (A  .GE.  O.DO)  ) 

+  GOTO  900 

Determine  temperature  and  orientation  of  detector  surface. 


N=1 

I F  (  -ZDN  .LT.  .ID-10  )  GOTO  220 
N=2 

IF  (  ZDN+ZCC  .GT.  .ID-10  )  GOTO  220 
N=3 

IF  (  DABS  (ZDN+ZCC)  .LT.  .10-10  )  GOTO  220 
N=4 

IF  (  ZDN+ZCCB  .GT.  .10-10  )  GOTO  220 
N=5 

CALL  TOCCS (N.PHIDN.TDS.ODS) 

Find  the  end  point  of  optical  path. 

IF  (  THEN  .GT.  .ID-10  )  GOTO  320 
IF  (  RHODN  . LE.  RCB  )  GOTO  310 
N=3 

GOTO  330 
N=5 

GOTO  340 

IF  (  DABS  (THEN-PI)  .GT.  .ID-10  )  GOTO  400 
N=1 

I F  (  RHODN  .LT.  RCB  )  GOTO  3^0 

REOP=ZCC 

GOTO  800  • 

REOP=ZCCB 
GOTO  800 

IF  (  RHODN  .GT.  .10-10  )  GOTO  410 

R1 =RP/STHEN 

R21=RCB/STHEN 

R22=0 . DO 

GOTO  450 

B1=A*A+RP*RP/RH0DN/RH0DN-1 .do 
B2=A*A+RCB*RCB/RH0DN/RH0DN-1 .DO 
R1=RHODN* (-A+DSQRT  (B1) ) /STHEN 
IF  (  B2  .GE.  O.DO  )  GOTO  420 
R2 1 =0  .  DO 
R22=0 . DO 
GOTO  450 

IF  (  RHODN  .GT.  RCB  )  GOTO  430 
R22=0 . DO 
GOTO  440 

R22-RH0DN*  (-A-DSQRT  (B2) ) /STHEN 
R21=RHODN*  (-A+DSQRT  (B2) ) /STHEN 
Z1  =-R1  -‘CTHEN+ZDN 
Z21=-R21*CTHEN+ZDN 
Z22=-R22*CTHEN+ZDN 

IF  (  -ZON-ZCC  .GT.  .10-10  )  GOTO  600 
IF  (  -Z 1  .GT.  O.DO  )  GOTO  520 
REOP=-ZON/CTHEN 
N-1 
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GOTO  800 

520  IF  (  -Z1  .GT.  ZCC  )  GOTO  5^0 
530  RE0P-R1 
N=2 

GOTO  700 

540  IF  (  (-Z21  .GT.  ZCC)  .AND.  (-Z22  .LT.  ZCC)  )  GOTO  550 
REOP-  (ZCC+ZDN) /CTHEN 
N=3 

GOTO  800 

550  IF  (  -Z21  .GT.  ZCCB  )  GOTO  580 
570  RE0P-R21 
N=4 

GOTO  700 

580  REOP- (ZCCB+ZDN) /CTHEN 
N=5 

GOTO  800 

600  IF  (  -Z21  .GT.  ZCCB)  GOTO  580 

620  IF  (  -Z21  .GT.  ZCC  )  GOTO  570 

61*0  IF  (  -Z1  .LT.  O.DO  )  GOTO  510 

660  GOTO  530 

700  IF  (  DABS  (RHODN)  .GT.  .10-10  )  GOTO  750 

720  PHIN— ZETN 

GOTO  800 

750  XEOP-REOP*STHEN/RHODN 

SPDN=DS I N  (PDN) -XEOP*SZETN 
CPDN-DCOS  (PDN) +XEOP*CZETN 
PH  I N-DATAN2  (SPDN.CPDN) 

800  CALL  TOCCS  (N.PHIN.TES.OES) 

REOP-DABS  (REOP) 

900  RETURN 
END 


C»> - «< 

C>»  «< 

C»>  Orientation  of  Combustion  Chamber  Surface  <« 

C»>  «< 

C»> - «< 


SUBROUTINE  TOCCS  (N , PH  I N ,TS ,0S) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKRD/TSCH ,TSCW,TSP,TSB , ES 
COMMON  /BLKRE/THEN , ZETN , STHEN , CTHEN , SZETN , CZETN 
IF  (  N-2  )  100,200,300 
100  TS-TSCH 
GOTO  700 
200  TS-TSCW 
GOTO  800 

300  IF  (  N-4  )  400,500,600 
1*00  TS=TSP 

GOTO  700 
500  TS-TSB 

GOTO  800 
600  TS-TSB 
700  OS-DABS  (CTHEN) 

RETURN 

800  OS-DABS  (STHEN*DCOS  (ZETN+PH I N) ) 

RETURN 
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END 


C»> - <<< 

C»>  «< 

C»>  RHF  along  an  OPTICAL  PATH  «< 

C»>  «< 

c>>> - «< 


DOUBLE  PRECISION  FUNCTION  RHFR (NOP , PDN , REOP) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR1/PI  ,PIST,PH* 

COMMON  /BLKR2/CFRE  (40)  , ERF C , ERFCL 

COMMON  /BLKR3/FVF,FVA,FVB.FVC,FAF,FAA,FAB,FAC,XHF,XHA,XHB,XHC, 
+  XCF , XCA , XCB , XCC 

COMMON  /BLKR7/REFL,REFTP,REFQ,REFW,REFWL,CVLWL 
COMMON  /BLKR8/OL,P,PE,T,TO,W 
COMMON  /BLKRA/RHOD , PH  I D , ZD , R I , STHE I N , CTHE I N 
COMMON  /BLKRB/THE , ZET 
COMMON  /BLKRC/TA.TB.TC 
COMMON  /BLKRD/TSCH,TSCW,TSP,TSB,ES 

COMMON  /BLKRF/QRS  (10)  .RHFRNR.QRSNR  (10)  .RHFRP.QRSP  (10) 

DIMENSION  CKO  (1 0)  ,Q0S(10)  ,QES(10)  ,  QRY  (5)  ,TAU  (5)  ,TAUR(10) 
DIMENSION  RS  (5) .PHIS  (5) 

RHFR-O.DO 

RHFRP-O.DO 

RHFRNR=O.DO 

NSUB-0 

DO  100  1=1,10 
QRS  ( I ) =0 . DO 
QRSNR  (I) =0.00 
100  QRSP  ( I } =0 .DO 
C»> 

C»>  Find  range  of  optical  path  and  surface  temperature. 

C»> 

CALL  ROP (REOP.TDS.TES.ODS.OES.NEOP) 

IF  (  REOP  .GT.  O.DO  )  GOTO  120 
RETURN 

120  DX=0 . 2D0 
DWN=0.7D0 
R1 =0 . DO 
R2-RE0P 

C>» 

C»>  Compute  surface  emission  power 
C»> 

TESl»=TES*TES*TES*TES 

TDS1*=TDS*TDS*TDS*TDS 

A-O.DO 

YE1-0.D0 

YD1-0.D0 

EE1-PI4/.15D2 

ED1-PI4/.15D2 

DO  150  1=1,9 

IF  (  TES  .GT.  .10-10  )  GOTO  1 30 
QES  (l)-O.DO 
GOTO  150 
YE1-YE1+DWN/TES 

IF  (  TDS  .GT.  .ID-10  )  GOTO  UO 


130 


QDS  ( I ) -0  -  DO 
GOTO  150 

140  YD1-YD1+DWN/TDS 

EE2-RIF  (3.0.D0.YE1) 

ED2-RI  F  (3.0.D0.YD1) 

QES ( I ) =TES4* . 15D2/P 1 4*  (EE  1 -EE 2) *OES 
QOS ( I ) =TDS4* . 1 5D2/P I  4* (EDI -ED2) *ODS 
EE1-EE2 
150  ED1=ED2 

QES (10)  *TES4*. 15D2/P I 4ftEE2*OES 
QOS  (10) *TDS4ft . 1 5D2/P I 4*ED2*OOS 
160  DO  170  1-1,10 
CKO  ( I ) =0 . DO 

IF  (  I  .GT.  5  )  GOTO  170 
QRY  (|)=0.D0 
TAU  ( 1) -0 . DO 
170  TAUR  (I) *0.00 
R 1*0. DO 

CALL  CTDL (THE , ZET) 

STHE-DS I N  (THE) 

IF  (  RHOD  .LE.  O.DO  )  GOTO  180 
NSUB*REOP*STHE/RHOD/DX 
180  IF  (  NSUB  .LT.  NOP  )  NSUB=NOP 
IF  (  NSUB  .LE.  50  )  GOTO  200 
NSUB-50 

R2-DX*RH0D*NSUB/STHE 
200  DR*  (R2-R1) /NSUB 
RI=R1+DR*.5D0 
OL-DR*REFL*CVLWL/REFWL 

C>» 

C»>  Heat  flux  integration  along  optical  path. 

C»> 

DO  500  1*1 .NSUB 

CALL  CTDL  (THE, ZET) 

210  SPH I D=DS I N  (PH  I D) 

CPH I D=DCOS  (PH ID) 

IF  (  DABS  (RHOD)  .LT.  .10-10  )  RHOO=O.DO 
IF  (  DABS  (ZD)  .LT.  .10-10  )  ZD=O.DO 

C>» 

C»>  Coordinate  transformations 
C»> 

CALL  CT (FVF , FVA, FVB, FVC, FVDL.RCFV.RWFV) 

CALL  CT (F AF , F AA, F AB , FAC , F ADL , RCF A ,RWF A) 

CALL  CT (XHF , XHA , XHB , XHC , XHDL , RCXH.RWXH) 

CALL  CT (XCF.XCA.XCB.XCC.XCDL, RCXC ,RWXC) 

C»> 

C»>  Compute  combined  volume  absorptance. 

C»> 

ZFA- (R! -RCFA) /RWFA 
ZFA-FADL-ZFAftZFA 

IF  (  DABS  (ZF A)  .GT.  .  1  72703  )  ZFA-ZF A/DABS  (ZFA) * . 1  7 2 7 D 3 
FA-DEXP(ZFA) 

IF  (  FA  .LT.  . 2D0  )  GOTO  500 
T-AFT(FA) 

ZFV* (R I -RCFV) /RWFV 


ZXH- (R I -RCXH) /RWXH 
ZXC-(RI-RCXC) /RWXC 
ZXH-XHDL-ZXH*ZXH 
ZXC=XCDL-ZXC*ZXC 
ZFV=FVDL-ZFV*ZFV 

IF  (  DABS  (ZFV)  .GT.  .1727D3  )  ZFV=ZFV/DABS  (ZFV) *. 1 727D3 
FV=0 . DO 

IF  (  DABS(FVF)  .GT.  O.DO  )  FV=DEXP  (ZFV) 

IF  (  DABS(ZXH)  .GT.  .172703)  ZXH=ZXH/DABS  (ZXH) * . 1 727D3 
XH=O.DO 

IF  (  DABS(XHF)  .GT.  O.DO  )  XH-DEXP  (ZXH) 

IF  (  DABS(ZXC)  .GT.  .172703  )  ZXC=ZXC/DABS  (ZXC) * . 1 727D3 
XOO.DO 

IF  (  DABS  (XCF)  .GT.  O.DO  )  XC-DEXP (ZXC) 

CALL  CVA  (FV,XH,XC,CKO,FVDL,XHDL,XCDL) 

.  IF  (  FVF  .NE.  O.DO  )  GOTO  212 
RWFV-RWXH 
RCFV-RCXH 

212  Z2*(RI+DR*.5D0-RCFV)/RWFV 

Z1« (Rl -DR*.5DO-RCFV) /RWFV 
214  WNI-O.DO 

ERFZ1-ERRF  (Z1) 

ERFCZ1-ERFC 
ERFLZ1-ERFCL 
ERFZ2-ERRF  (Z2) 

ERFCZ2-ERFC 

ERFLZ2-ERFCL 

C»> 

C»>  integrate  along  optical  path  at  each  spectra-1  band. 

C>» 

DO  400  J-1,10 

IF  (  (DABS  (FVDL)  .GT.  .1727D3)  -OR.  (DABS  (XHDL)  .GT. 

+  . 1 727D3)  -OR.  (DABS(XCDL)  .GT.  .1727D3)  )  GOTO  220 

TAUW-.5D0*PIST*CK0 (J)  *RWFV*REFL*CVLWL/REFWL 
TAUC-TAUR  (J) -TAUWAERFZ1 
TAURU-T  AUW*ERF  Z2+T  AUC 
IF  (  DABS  (ERFCZ1)  .GT.  1 .0-15  )  GOTO  230 
TAUC-TAUR (J) +TAUW*ERFCZ1 
TAURU-TAUC-TAUW*ERFCZ2 
GOTO  230 

220  TAUW1-CK0  (J) +ERFLZ1 

IF  (  DABS  (TAUW1 )  .GT.  .1727D3  )  TAUW1-DABS (TAUW1) / 

+  TAUW1*. 1727D3 

TAUW1-DEXP (TAUW1) *.500*PtST*RWFV*REFL*CVLWL/REFWL 
TAUW2-CK0(J)+ERFLZ2 

IF  (  DABS  (TAUW2)  .GT.  .  1727D3  )  TAUW2-DABS (TAUW2)  / 

+  TAUW2*.1727D3 

TAUW2-DEXP (TAUW2) * . 5D0*P I ST*RWFV*REF L*CVLWL/REFWL 
TAUC-TAUR  (J) +TAUW1 *DABS  (Z1 ) /Z1 
TAURU-TAUC-TAUW2*0ABS  (Z2) /Z2 
230  IF  (  TAURU  .LE.  TAUR  (J)  )  GOTO  400 

CALL  GNI X  (TAURU, TAUR (J) ,TAU) 

DO  300  K-1,5 

IF  (  J  .GT.  1  )  GOTO  260 

IF  (  (DABS  (FVDL)  .GT.  .  1727D3)  -OR.  (DABS  (XHDL)  .GT. 
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.172703)  .OR.  (DABS(XCDL)  .GT.  .1727D3)  )  GOTO  240 

TAUS-  (TAU  (K) -TAUC) /TAUW 

TAUSC-TAUS/OABS  (TAUS) -TAUS 

IF  (  DABS  (ERFCZ1)  .GT.  1.D-15  )  GOTO  250 

TAUD-TAUC-TAU  (K) 

.TAUSL-DLQG  (DABS  (TAUD) )  -DLOG  (TAUW) 

TAUSC-1 .0-75 

IF  (  DABS  (TAUSL)  .LT.  .  1727D3  )  TAUSC=DEXP  (TAUSL) 

TAUSC-TAUD/DABS (TAUD)  *TAUSC 

TAUS-TAUD/DABS (TAUD) *.99999900 

GOTO  250 

YY-TAUC-TAU  (K) 

TAUS-YY/DABS  (YY) *.9999999D0 
TAUSC=1 .0-75 

XX-.5D0*PIST*RWFV*REFL*CVLWL/REFWL 
YY-DABS (YY) 

IF  (  YY  .LE.  0 . DO  )  YY=1 .0-75 
TAUSL-DLOG  (YY) -CKO (J) -DLOG (XX) 

RS  (K)  -FRE  (TAUS, TAUS C .TAUSL)  *RWFV+RCFV 
IF  (  RWFA  .GT.  1  .070  )  RCFA-RS  (K) 

ZFA-  (RS  (K) -RCFA) /RWFA 
ZFA-FADL-ZFA*ZFA 

IF  (  DABS (ZFA)  .GT.  .1727D3  )  ZFA-ZFA/DABS (ZFA)  *.172703 
FA=DEXP (ZFA) 

IF  (  FA  .GT.  . 2D0  )  GOTO  280 
QRY  (K) *0 . DO 
GOTO  300 
T-AFT  (FA) 

A-T*TAU  (K) 

Y1-WN1/T 

Y2- (WN1+0WN) /T 

RIF2-0.D0 

IF  (  J  .LT.  10  )  RIF2-RIF  (L.A.Y2) 

R I F 1 =R I F  (A,A,Y1) 

QRY (K)  *T *T *T *T *T * (R I F 1 -R I F2) * . 1 502/P  I U 
CONTINUE 

QRS  (J)-QRS  (J)+GNI  (TAURU.TAUR  (J)  ,QRY) 

TAUR  (J) -TAURU 
WN1-WN1+DWN 
RI-RI+DR 

IF  (  R2  .GE.  REOP  )  GOTO  520 

R1-R2 

R2-RE0P 

NSUB-NOP 

GOTO  200 

Net  radiation  heat  flux  on  detection  surface,  (mean  reflection) 

RFL«1 .DO-ES 
DO  600  1-1,10 

X-.35DO+.7DO*(l-1) 

QRSP  ( I ) -QRS  ( I ) *Q0S 
RHFRP-RHFP.P+QRSP  (I) 

QRSNR  ( I ) -ES*OOS*QRS ( I ) -QDS ( I ) 

RHFRNR-RHFRNR+QRSNR ( I ) 


EXT-O.DO 

IF  (  TAUR(I)  .LT.  1.02)  EXT=DEXP (-TAUR  ( I ) *X) 

BES=(QES (l)+RFL*OES*QRS (I )+RFL*OES*EXT*  (QDS  (l)+RFL*ODS 
+  *QRS  (I))) / (1 . DO-RFL*RFL*OES*ODS*EXT*EXT) 

QRS  ( I ) =ES*OOS*  (QRS ( I )  +BES*EXT) -QDS ( I ) 

TAUR  ( I ) =EXT 
600  RHFR-RHFR+QRS  ( I ) 

RETURN 

END 


C>» - «< 

c»>  «< 
C»>  COMBINED  VOLUMETRIC  ABSORPTANCE  «< 
C»>  «< 
C»>  OL  -  dimensionless  optical  length  (ref:6u)  <« 
C»>  T  -  dimensionless  temperature  (ref:2400K)  «< 
C»>  P  -  dimensionless  pressure  (ref:1atm)  «< 
C»>  «< 
C»> - <<< 


SUBROUTINE  CVA (FV,XH20,XC02,CK0,FVDL,XHDL,XCDL) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR8/0L,P,PE,T,T0,W 
DIMENSION  CKO  (10)  ,VU  (11)  ,VL  (11) 

1  FORMAT  ('  ***  CVA  error  :  (fv,XH20,XC02)  =  ’ , 3E1 0 . 3) 

2  FORMAT  ( 1  ***  CVA  error  :  ( I -EMI S-CKO)  -  ' . 1 5.2E15-3) 

IF  (  (FV  . LE .  1 . DO)  .AND.  (XH20  .LE.  1.DO)  .AND. 

+  (XC02  .LE.  1  . DO)  )  GOTO  100 

WRITE  (6.1) FV.XH20.XC02 
STOP 

100  CALL*  GAS (XH20.XC02 jVU.VL) 

IF  (  VU(1)  .LT.  .10-10  )  VU(1)-0.D0 
J-1 

NWN-10 

DWN=7 . DO/NWN 
WN-O.DO 
DO  500  1=1, NWN 
A=0.D0 

150  IF  (  (I  .GE.  NWN)  .AND.  (J  .LE.  11)  )  GOTO  200 

IF  (  (VL  (J)  .GE.  WN+DWN)  .OR.  (J  .GT.  11)  )  GOTO  400 
200  Cl =WN 

IF  (  VL  (J)  .GT.  WN  )  Cl =VL  (J) 

C2-VU  (J) 

IF  (  VU  (J)  .LT.  WN+DWN  )  GOTO  300 
A=A+WN+DWN-C1 
VL  (J) -WN+DWN 
GOTO  400 
300  A-A+C2-C1 

J-J+1 
GOTO  150 

400  WNC-WN+ . 5D0*DWN 

WL-1 . DO/WNC 
EMI S-A/DWN 

IF  (  EMIS  .GE.  O.DO  )  GOTO  440 
420  WRITE (6,2) I , EMIS, CKO  (I ) 

STOP 

440  IF  {  EMIS  .GE.  1.D0-1.D-15  )  EMI S-1 .DO-1 .D-15 
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CKO { I ) -SOOTA (WL , F V) -DLOG (1 . DO-EMI S) /WNC/OL 

IF  (  CKO  (I)  .LT.  0 . DO  )  GOTO  1*20 

IF  (  CKO  (I )  .EQ.  0.00  )  GOTO  500 

IF  (  FV  .NE.  0 .DO  )  GOTO  1*50 

IF  (  DABS  (XHDL)  .LE.  .1727D3  )  GOTO 

CKO  { I ) =DLOG  (CKO  ( I ) ) +XHDL-DLOG  (XH20) 

GOTO  500 

1*1*5  CKO  (I)  =CKO  (I)  *DEXP  (XHDL)  /XH20 
GOTO  500 

450  IF  (  (DABS  (FVDL)  .LE.  .172703)  -AND.  (DABS  (XHDL)  .LE. 

+  .172703)  -AND.  (DABS(XCDL)  .LE.  .172703)  )  GOTO  460 

CKO ( I ) -DLOG (CKO (I) )+FVDL-DLOG (FV) 

GOTO  500 

1*60  CKO  (l)-CKO(l)  *DEXP  (FVDL)  /FV 

500  WN-WN+DWN 
RETURN 
END 

QAAAAAAAAA  A  AAA  A  A  it  it  it  A  A  A  A  A  A  A  A  it  i;  A  A A  A  *  A  A  A  A  A  A  A  A  A  A  i;  A  A  >’c  A  A  A  A  A  A  A  A  A  A  A  A  A  s':  A  A  A  A  A  A  A  A  A 


C*AA  AAA 

C***  SOOT  RADIATION  *** 

c***  AAA 

C*A*  AAA 

C AAA  AAA 


CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

c»> - - - 

c»> 

C»>  MEAN  SPECTRAL  SOOT  ABSORPTIVITY 

C»> 

C»>  .  1966  -  Howarth,  Forster  and.  Thring 

C>»  1969  “  Dalzell  and  Sarofim 

C»>  1981  -  Lee  and  Tien 

C»> 

C»>  WL  -  dimensionless  wave  length  (  ref:6u 

C»>  FV  -  soot  volume  fraction 

C»>  SOOTA  -  mean  soot  absorptanoe 

C»> 

c»> - 

DOUBLE  PRECISION  FUNCTION  SOOTA (WL.FV) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR1/PI .PIST.PI4 
IF  (  WL  .GE.  0.D0  )  GOTO  200 

100  WRITE(6,101) 

101  FORMAT  (IX, 'ERROR  in  SOOTA  function.') 

STOP 

200  CALL  SOC(WL.X.Y) 

SOOTA-l8.DO*PI*Y*FV/  ((X+2.D0) **2+Y*Y) 

RETURN 
END 


C»> - <<< 

c>»  «< 

C»>  OPTICAL  CONSTANTS  OF  SOOT  «< 

C>»  «< 

C»>  1981  -  Lee  and  Tien  «< 

C»>  «< 

C>»  -  dimensionless  operation  «< 


«< 

«< 

«< 

«<; 

«< 

«< 

«< 

)  «< 
«< 
«< 
«< 
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electronic  dispersion  function 

dispersion  constants,  F ,  G  £  WN  (Lee  £  Tien) 


SUBROUTINE  S0C(WL,X,Y) 

IMPLICIT  REALA8  (A-H.O-Z) 

COMMON  /BLKR1/PI ,P I  ST, PI  4 
COMMON  /BLKR6/F  (3)  ,G  (3)  ,WN  (3)  , VLGT 
X=1 .DO 
Y=0 .DO 

W=2.DO*PI*VLGT/WL 
DO  200  1=1,3 

A=WN  ( I ) *WN ( I ) -W*W 
B=A*A+W*W*G ( I ) *G ( I ) 

IF  (  DABS  (B)  .LT.  1.D-75  )  GOTO  900 
100  X=X+F  ( I ) *A/B 

Y=Y+F  (I)  *G  (I)  *W/B 
200  CONTINUE 
900  RETURN 
END 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 
CAAA  AAA 

C***  GAS  RADIATION  aaa 

CAAA  AAA 

C***  1972  -  Edwards  and  Balakrishnan  *** 

C***  1983  "  Chang  and  Rhee  aaa 

C**A  aaa 

C aaa  X???  -  molar  fraction  of  gas  ???  *** 

C aaa  VU  -  dimensionless  upper  bound  of  gas  bands  *** 

C aaa  VL  -  dimensionless  lower  bound  of  gas  bands  *** 

C*a*  AAA 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

SUBROUTINE  GAS  (XH20, XC02 , VU, VL) 

IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  AA1  (5)  ,AA2  (6)  ,VU  (11)  ,VL  (11) 

1  FORMAT  ('  XH20  =  1 , El  2.4,'  XC02  =  ’.E12.U) 

I F  (  (  XH20  .GE.  0 .00  )  .OR.  (  XC02  .GE.  0.D0  )  )  GOTO  100 

WRITE  (6,1) XH20.XC02 

STOP 

100  VU  (2) =. 160D4 
VU  (3)=.37&D4 
VU  (4) =.53504 
VU  (5) -.72504 
VU  (6) =.66703 
vu  (7) -.960D3 
VU  (8)  -  .10604 
VU(10)=. 36604 
VU(11)=. 52004 
CALL  AH20  (XH20, AA1 ) 

CALL  AC02  (XC02.AA2) 

VU(1)  =  A  A 1  (1)  A. 500 

VU  (9) = . 241 0D4-AA2 (4)  A.5D0 

DO  200  1=1,5 

VU ( I )  =VU ( I ) +AA1  ( I ) * .500 
200  VL ( I )  =VU ( I )  -AA1 ( I ) 


DO  300  1=1,6 

VU ( I  +5)  =VU ( I +5) +AA2  ( I ) * . 50 o 
300  VL  (1+5)  *VU  (1+5)  -AA2  (I) 

DO  1*00  1=1,10 
11=1+1 

DO  400  J=I1 ,11 

IF  (  VL  (I)  .LE.  VL{J)  )  GOTO  400 
TEMP=VL  (I) 

TEMPI =VU  ( I ) 

VL  (I)  =VL  (J) 

VU.(I)=VU  (J) 

VL  (J)  =TEMP 
VU (J)  =TEMP1 
400  CONTINUE 

DO  500  1=1,10 

I  F  (  VU(I)  .LE.  VL  (1+1)  )  GOTO  500 
I F  (  VU  ( I )  .LT.  VU  (1+1)  )  GOTO  420 
VU  ( I +1 ) =VU  ( I ) 

420  VU  (I)  -VL  (1+1) 

500  CONTINUE 

DO  600  1=1,11 

VU(l)=VU(l)*.6D-3 
600  VL  (I)  =VL  (I)  *.6D-3 
RETURN 
END 


C>» - «< 

c>»  «< 

C»>  BAND  ABSORTANCE  of  WATER  «< 

C»>  «< 

C»>  1973  “  Edwards  and  Balakrishnan  «< 

C»>  «< 

C»> - <<< 


SUBROUTINE  AH20(X,AA) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  AA(5),B(3),ND(3)  ,NG(3)  ,V(3) 

COMMON  /BLKR8/0L,P,PE,T,T0,W 

DATA  NG/1 ,1,1/, V/3.652D3.1 .59503-. 3- 75603/ 

I F (  (X  ,GT.  0 . DO)  .AND.  (T  .GT.  0.D0)  )  GOTO  200 
DO  100  1=1,5 

100  AA ( I ) =0 . DO 
RETURN 

200  DST=0SQRT  (T/TO) 

PE=P+P*X*(.86D1/DST-.5D0) 

W=P*X*OL/T/1 .824D3 
B  (1)=.52D4*DEXP  (-.176D2/DST) 

B  (2)  =  . 1431100 
B(3)=.284D2 
ND  (1 )  =0 
ND  (2)  =0 
ND(3)-0 

CALL  A1  (V,NG, B,ND, AA (1) ) 

B  (1)  =  .4l  2D2 
B  (2) -.91*270-1 
B  (3) *-56402 
ND (2) -1 
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CALL  A1  (V, NG , B, ND , AA  (2) ) 

B  (1)  = .  1 900 
B  (2)  =.  1 321 9D0 
B  (3)  =  .602 
ND  (2)  =2 

CALL  PSIPHI  (TO.V.NG.ND.PSIO.PHIO) 

CALL  PSIPHI  (T.V.NG, NO, PS  I , PHI) 

T1-B  (1)*PSI/PSI0 
TA=B  (2)  *PH I /PH  I  0 
B  (1)-.23D1 
ND  (1)=1 
ND  (2)  *0 

CALL  PSIPHI  (TO.V.NG.ND.PSIO.PHIO) 

CALL  PSIPHI  (T,V,NG,ND,PSI , PHI) 

T2-B  (1)  *PS  I /PS  I  0 
T5*B (2)  ftPHI/PHI 0 
B  (1)  -.224D2 
ND  (1 )  =0 
ND(3)=1 

CALL  PSIPHI  (TO.V.NG.ND.PSIO.PHIO) 

CALL  PSIPHI  (T.V.NG. ND.PSI .PHI) 

T3“B  (1) *PS I /PS  I  0 
T6-B (2)  *PH I /PH  I  0 
C1-T1+T2+T3 
C3*B  (3) *DST 

C2=2 .DO*DSQRT (C3) * (DSQRT (T1*TL) +DSQRT (T2*T5) +DSQRT (T3*T6) ) 
CALL  BANDA  (Cl .C2.C3.AA  (3) ) 

B(1)-3.D0 
B(2)».8l69D-1 
B(3)“.U31D2 
ND (2) *1 

CALL  A1  (V.NG.B.ND, AA  {k) ) 

B(1)*.25D1 
B(2)=.11628D0 
B  (3)  32D2 

ND  (1)=1 
ND  (2)  =0 

CALL  A1  (V, NG , B , ND , AA  (5)) 

RETURN 
EN 


PS  I  and  PHI  functions  of  GAS  EMISSION 


1968  -  Weiner  and  Edwards 


SUBROUTINE  PS  I  PH  I  (T, V, NG , NO , PS  I , PH  I ) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  V  (3)  .NG  (3)  ,U(3)  ,N0(3)  ,NV(3) 

COMMON  /BLKR7/REFL, REFTP, REFQ.REFW.REFWL.CVLWL 
IF  (  (NO  (1 )  .NE.  0)  .OR.  (NO  (2)  .NE.  0) 

+  .0R.  (ND  (3)  -NE.  0)  )  GOTO  100 
PSI-1 .00 

PHI-1 .DO/OSQRT  (T) 


RETURN 

M=3 

UD=0 . DO 
XD=0 . DO 
XY=0 . DO 
YN=0 .  DO 
DO  400  K=1 ,M 

U  (K)  =V(K)  *REFWL*1  .D-4/T 
NV  (K)  =0 

IF  (  NO  (K)  .LT.  0  )  NV  (K)  =  I  ABS  (N  D  (K) ) 

NDK= I ABS  (ND  (K) ) 

Y1=0 . DO 

X=DEXP  (-U  (K)  *NV(K))  /PER  (NG  (K)  -1  ,  NG  (K)  -1) 
X1=X*PER  (NV  (K)  +NG  (K)+NDK-1  , NG  (K)+NDK-1) 
X2=X*PER (NV (K)  +NG  (K) -1 , NG  (K) -1) 

I F (  XI  .LT.  1 . D  —  1 0  )  GOTO  300 

XY=XY+X1 

XD=XD+X2 

Y1=Y1+DSQRT  (XI) 

NV(K)=NV(K)+1 

IF  (  NV  (K)  .LT.  100  )  GOTO  200 

YN=YN+Y1*Y1 

UD=UD-U  (K) *ND  (K) 

CONTINUE 

PS  1 =  (1 .DO-OEXP  (UD) )  *XY/XD 
PH  I =YN/XY/DSQRT (T) 

RETURN 

END 


APPROXIMATION  to  GAS  BAND  ABSORPTANCE 


1964  "  Edwards  and  Menard 


SUBROUTINE  BANDA  (Cl , C2 , C3 , A) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR8/OL,P,PE.T,TO,W 

FORMAT (T20, ******  ERROR  IN  BANDA  ****** ,/,T20,5E12. 4) 

I F  (  (Cl  .LT.  0 . DO)  .OR.  (C2  .LT.  0.00)  .OR.  (C3  .LT.  O.DO)  ) 
+  GOTO  200 
IF  (  W*PE  )  200,100,300 
A=0 . DO 
RETURN 

WRITE  (6,1) Cl ,C2,C3,PE,W 
STOP 

E-C2*C2*PE/ (4.DO*C1*C3) 

I F  (  E  .GT.  1 .DO  )  GOTO  400 

TEST1=E*C3 

TEST2*C3* (2 . DO— E) 

A=C1 *W 

IF  (  A  .LE.  TEST  1  )  GOTO  900 
A=C2*DSQRT (W*PE) -E*C3 
IF  (  A  .LE.  TEST2  )  GOTO  900 
X=E*C1*W/C3 


A=C3*(DL0G  (X)  +2.D0-E) 

GOTO  900 
400  TEST=C3 
A=C1*W 

I F  (  A  . LE .  TEST  )  GOTO  900 
X=C1*W/C3 

A=C3* (CLOG  (X) +1 .DO) 

900  RETURN 
END 


C»> - <<< 

c>>>  <<< 

C»>  GAS  BAND  ABSORPTANCE  (regular  bands)  <« 

C»>  «< 

c»> - <<< 


SUBROUTINE  A1  (V, NG , B , ND , A) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  B  (3)  ,ND  (3)  ,V(3)  ,NG  (3) 
COMMON  /BLKR8/OL,P,PE,T,TO,W 
CALL  PSIPHI  (TO,V,NG,ND,PSIO,PHIO) 
CALL  PSIPHI  (T, V, NG , ND , PS  I , PH  I ) 

C3=B  (3)  *DSQRT  (T/TO) 

C1=B  (1)  *PS  I /PS  I  0 

C2=2 . DO*DSQRT  (Cl *C3*B  (2) *PH I /PH  I  0) 
CALL  BANDA  (Cl ,C2,C3, A) 

RETURN 

END 


C»> - <<< 

C»>  <<< 

C»>  ;  BAND  ABSORPTANCES  of  CARBON  DIOXIDE  •  «< 

C»>  <« 

C»>  1973  “  Edwards  and  Balakrishnan  «< 

C»>  «< 

C»> - «< 


SUBROUTINE  AC02(X,AA) 

.  IMPLICIT  REAL*8  (A-H,0-Z) 

DIMENSION  AA  (6)  ,B  (3)  .NO  (3)  . NG  (3)  ,V(3) 
COMMON  /BLKR8/0L,P,PE,T,T0,W 
DATA  NG/ 1 ,2,1/, V/ 1 . 35 1 03 , 6 . 6702 , 2 . 39&D3/ 
DO  100  1=1,6 
100  AA  ( I ) =0 . DO 

W=P*X*0L/T/7. 46202 

IF  (  W  .LE.  0 . DO  )  RETURN 

PE=  (P+P*X*0 . 3D0)  **0 . 7D0 

B  (1 )  = .  1 9D2 

B (2) *.61570-1 

B (3) = - 1 27D2 

ND  (1)  =0 

ND (2) =1 

ND  (3)  =0 

CALL  A1  (V, NG , B , ND , AA (1) ) 

PE=  (P+P*X*0.3D0) **0 . 8D0 
B  (1)  =0.3600 
B (2) =.40170-1 
B(3)  =  - 13^02 
NO  (1)=-1 
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ND  (2)  *0 
ND(3)“1 

CALL  A1  (V, NG , B , ND , AA (2)) 

B(1)=.24D0 
B (2) = . 1 1 888D0 
B  (3)  =.10102 

CALL  A1  (V, NG , B , ND , AA (3) ) 

B (1) =. 1 1 D3 
B (2)  = . 24723D0 
B  (3) *-11202 
ND  (l)  =0 

CALL  A1  (V,NG,B,ND, AA (4) ) 

PE=  (P+P*X*0 . 3D0)  **0 . 6500 
B  (1 )  =1* .  DO 
B (2) =0.1 33^1 

B (3) =23.5 

ND  (1)=1 

CALL  A1  (V,NG,B,ND,AA(5)) 

B  (1) =.66D-1 
B  (2) =.3930500 
B (3) *.3^502 
ND  (1)=2 

CALL  A1 (V, NG , B , ND , AA (6)) 

RETURN 

END 

c  it  it  it  it  it  )'c  *  *  it  it  it  *  it it  * it  i;  it  *  i;  it  it  it  it  it  it  it  it  it  it  it  it it it 


it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it  it 


it  it  it  it  it  it  it  it 


C  it  it  it  it  it  it 

Ci<*i<  Coordinate  Transformation  ’  *** 

Qititit  >  ititlt 

C***  Fuel -plume  < - >  Optical  path  *** 

C  ***  it  it  it 


C*********ft**********************ft*******rtft***:ft*rtrtrt*ft***y:*rtrtrt*ft*ft*ft*ftrt* 

SUBROUTINE  CT (FF . A.B.C.FOL.RC.RW) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR4/RHOON , PH  I  ON . ZDN , THE  I N ,OMES , RNZ , VS , NF P 
COMMON  /BLKRA/RHOD . PH  I D , ZD . R I , STHE I N , CTHE  I  N 
COMMON  /BLKRB/THE , ZET 

1  FORMAT  (3X » ' ***  ERROR  in  CT  -  (THE : ZET:RI:A1:B1:C1)  =  ',/,27X. 

+  6E12.3) 

STHE=DS i N  (THE) 

CTHE=DCOS  (THE) 

SZET=DSIN  (ZET) 

CZET=DCOS  (ZET) 

SPH I D=DS I N  (PH  I D) 

CPHID-DCOS(PHID) 

IF  (  FF  .EQ.  0.00  )  GOTO  20 

IF  (  (A  .NE.  0 . DO)  .OR.  (B  .NE.  O.DO)  .OR.  (C  .NE.  O.DO)  ) 

+  GOTO  40 

20  FDL—1.00 

IF  (  DABS  (FF)  .GT.  O.DO  )  FDL=DLOG  (FF) 

RC=R  I 
RW=1 .070 
RETURN 

ALP-.850-1/RNZ 

PH  I  Si*  (OMES*CTHE I N/VS) * (1 . DO+ALP* .500) 


40 


40 


IF  (  OABS(RHOD)  .GT.  .ID-10  )  GOTO  100 


Radial  component  (RHO  term) 


SR=1 .DO 
X  I  =R  I 

RHO£*R I *STHE 
RHOA*RHOE 
RHOB-STHE 
RHOC*O.DO 
RHOO-O.DO 
RH01-STHE 
RHO2-O.D0 
RHO-RHOE 
GOTO  230 


I  I .  RHOD  <>  0 


SR-STHE/RHOD 
X I  -R I *SR 

ARHO=DCOS (PHID+ZET) 

IF  (  XI  .IT.  1.D1  )  GOTO  120 

RHOO-RHOD*ARHO 

RH01=RH0D 

RHO2-0.D0 

GOTO  200 

FA=DSQRT (1 . DO+2 . DO*ARHO*X l+X I *X I ) 

RHOA=RHOD*FA 
RHOB=RHOD* (ARHO+X I )  /F A 

RHOORHOD*  (FA*F A-  (ARHO+X  I )  *  (ARHO+X  I ) )  /FA/FA/FA 
RHOO=RHOA-RHOB*X I +RHOC*X I *X I  * . 5D0 
RH01*RH0B-RH0C*X I 
RH02=RHOC*.5D0 
RHO*RHOO+RHQ1 *X I +RH02*X I *X  I 


Angular  component  (PHI  term) 


IF  (  THE  .NE.  0.D0  )  GOTO  220 

RHOE=RHOD 

GOTO  230 


EXACT  TRANSFORMATION 


RHOE*RHOO*DSQRT (1 . DO+2 . DO*ARHO*X I +X I *X I ) 

X1-CPHID+XI*CZET 

X2»SPH I D-X I *SZET 


V  V  r*  / 


IF  (  (XI  .EQ.  0 . DO)  .AND.  (X2 
PH  I E=DATAN2  (X2.X1) 

ZE=ZD-R I *CTHE 


0. DO)  )  XI =1  -  DO 


,  ID-10  )  GOTO  290 


Approximated  Transformation  -  swirl  effect 


PH  I SA=PH I  SI *RHOE*  (1 . DO+ALP*RHOE*.5DO) / (1 . DO+ALP* -5D0) 

PH  I SB=PH I  SI*  (1 . DO+ALP*RHOE) / (1 . DO+ALP* . 5D0) *RHOB 

PH  I  SOPH  1  SI  *  A  L  P  /  (1  .DO+ALP*. 5D0)  *RHOB*RHOB+PH  I  SB*RHOC 

IF  (  (RHOD  .GT.  .ID-10)  .AND.  (  XI  .LT.  1.D1)  )  GOTO  250 

PH  I A=-ZET 

PH  I B=0 .DO 

PH  I OO .DO 

GOTO  300 

X3=DSIN  (PHID+ZET) 

IF  (  DABS (X3)  -GT.  .10-10  )  GOTO  280 
PH  I B=0 .DO 
PH  I C=0 . DO 

IF  (  RHO  .LT.  0 . DO  )  GOTO  260 
PHI A-PHID 
GOTO  300 
PH  I A=-PH I D 
RHO=-RHO 
GOTO  300 
PH  I A=PH I E 

IF  (  DABS  (XI)  .GT.  .ID-10  )  GOTO  290 

PHIB—CZET*CZET/X3 

PHIC=O.DO 

GOTO  300 

F  A=X2/X1 

FB—X3/X1/X1 

FC— 2.D0*FB*CZET/X1 

G=1 . DO+FA*FA 

PH  I B=FB/G 

PH  I C=  (-2 . DO*F  A*FB*FB+G*F  C)  /G/G 

PS  I A=PHI A-PHISA 

PS  I B=PH I B-PH I  SB 

PS  I C=PH I C-PH I  SC 

BPH I A*PS I A*PS I  A 

BPH I B=2 .DO*PS I A*PS I B 

BPH  I  02 .  DO*PS  I  A*PS  I  C+2 .  DO*PS  I  B*PS  I  B 

BPH I O-BPH I A-BPH I 8*X I +BPH I C*X I *X  I  * . 5D0 

BPH I 1»BPH I B-BPH I C*X I 

BPHI2-BPHIC*.5D0 


Coefficients  -  Approximated  Transformation 


A1«A*RH00+B*BPHI 0+C*ZD*ZD 
B1—  (A*RH01+B*BPH1 1)  *SR+2.D0*C*ZD*CTHE 
C1«  (A*RH02+B*BPH  I  2)  *SR*SR+C*CTHE*CTHE 
IF  (  Cl  .GT.  O.DO  )  GOTO  420 


RWgSac5qw 


42 


COMMON  /BLKRE/THEN , ZETN , STHEN , CTHEN , SZETN , CZETN 
ALP=.85D-1/RNZ 
X I =R I *STHEN 

RH01 =X I *X 1+2 .DOfcDCOS (PH I DN-ZETN)  *X I *RHODN+RHODN*RHODN 
RHO1=0SQRT  (RHOI) 

Y»RHODN*DS I N (PH  I  ON) -X I *DS I N (ZETN) 

X- RHODN*DCOS (PH  I DN) +X I *OCOS (ZETN) 

PHI 1=-ZETN 

IF  (  (DABS  (X)  .GT.  .ID-10)  .OR.  (DABS  (Y)  .GT.  .ID-10)  ) 

+PH 1 1=DATAN2  (Y,X) 

Z1=ZDN-R I *CTHEN 

JFP-DABS (PH  1 1 ) *NFP/2 .DO/P 

PH  I  FP=P 1*2. DO*JFP/NFP 

IF  (  PH  1 1  .LT.  O.DO  )  PH  I FP=-PH  I  FP 

ZETFP=-PH I FP 

CALL  CTI FP(ZETFP,THEIN, THEN, ZETN, THE, ZET) 

PDN=PH I DN-PH I FP 

XI- RHODN*DCOS (PDN) 

Y1=RHODN*DS I N (PDN) 

X2*X1 *CTHE I N-ZDN*STHE I N 
Y2-Y1 

ZD=X1 *STHE I N+ZON*CTHE I N 
RHOD-DSQRT (X2*X2+Y2*Y2) 

PH  I D=PDN 

IF  (  (DABS  (Y2)  .GT.  .10-10)  .OR.  (DABS  (X2)  .GT.  .ID-10)  ) 

+  PHID=DATAN2 (Y2.X2) 

RETURN 

END 

C***ft**rt***rt********rt*rt3V**rtrt3*rt*ft**#a*****rtft**fcyc*rt*********rt***5l:ftftrt*Yc**rt 
C**A  *** 

C*A*  2-D  TRANSIENT  CONDUCTION  *** 

C***  -  Finite  difference  approximation  in  space  *** 

C*a*  and  successive  d i f ferenc i at  ion  in  time  *** 


C*A*  2-D  TRANSIENT  CONDUCTION  A** 

C**A  -  Finite  difference  approximation  in  space  *** 

C*aa  and  successive  d i f ferenc i at  ion  in  time  *** 

CAAA  AAA 

C*A*  (CONFU . FOR)  AAA 

CAAA  AAA 

CAAA  ■kiti; 

C  is  kit  AAA 

Qkkk  AAA 

C**Yt  A  A  A 

CAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 

SUBROUTINE  MATRIC 
IMPLICIT  REAL*8  (A-H.o-Z) 

COMMON  /BLKC1/XC  (6 0)  ,YC(60)  ,C(60)  ,CK(60) 

COMMON  /BLKC2/A  (60,60)  ,B  (60,60) 

COMMON  /BLKC3/D  (60) , UBD  (60) , UO  (60) , UB (60) 

COMMON  /BLKC4/REFT,REF0T,REFC,REFK,REFLC,NPT,NU,NBC,NBCT,NLP 
COMMON  /BLKC5/L1  (80)  ,  L2  (80)  ,  L3  (80)  ,  1 1  (80)  ,  I  2  (80)  ,  I  3  (80)  , 

1  J1  (80)  ,J2(80)  ,J3(8o) 

COMMON  /BLKC6/REFQC,REFA,REFTK,REFU 


WRITING  FORMATS 


fl 

»*■ 


s 

m I 


ws 


SM 


s 

1 


FOL—  .1727D3 
RC-O.DO 
RW«1 .00 
RETURN 

IF  (  Cl  .LT.  .ID-10  ) 

X»A1-B1*B1*. 2500/Cl 

FDL-DLOG  (FF) -X 

RC*B1/C1*.5D0 

RW*1 . DO/DSQRT (Cl) 

RETURN 

END 


C1=. 10-10 


Spherical  Coordinate  Transformation 


SUBROUTINE  CT I FP  (ZETFP .THE  I N , THEO, ZETO.THE , ZET) 

IMPLICIT  REAL*8  (A-H.O-Z) 

DIMENSION  EX  (3)  ,  EY  (3)  ,EZ(3)  ,PT(3) 

STHE IN-OS  IN (THE  IN) 

CTHEIN-OCOS (THE  IN) 

CZETFP-OCOS (ZETFP) 

SZETFP-OS IN (ZETFP) 

EX  (1) “CTHE I N*CZETFP 
EX  (2) *CTHE I N*SZETFP 
EX (3) “STHE I N 
EY  (1)=-SZETFP 
EY  (2) “CZETFP 
EY(3)=O.DO 

EZ  (1)  “-STHE  IN*CZETFP 
EZ  (2) “-STHE I N*SZETFP 
EZ  (3) “CTHE I N 

PT  (1)  =DS  I  N  (THEO)  *OCOS  (ZETO)  * 

PT (2)  =DS I N (THEO) *DS I N  (ZETO) 

PT(3)  -DCOS  (THEO) 

CTHE-PT  (1)  *EZ  (1)  +PT  (2)  *EZ  (2)  +PT  (3)  *EZ  (3) 

STHE=DSQRT  (1 . DO-CTHE*CTHE) 

THE“DATAN2 (STHE.CTHE) 

X=PT  (1)  *EX  (1 )  +PT  (2)  *  EX  (2)  +PT  (3)  *EX  (3) 

Y-PT  (1 )  *EY  (1 )  +PT  (2)  *EY  (2)  +PT  (3)  *EY  (3) 

ZET-ZETO 

IF  (  (DABS  (Y)  .GT.  .10-10)  .OR.  (DABS  (X)  .GT.  .10-10)  ) 
+  ZET-DATAN2  (Y.X) 

RETURN 

END 


Coordinate  Transformation  of  Detector's  Location 


SUBROUTINE  CTDL  (THE ,ZET) 

IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON  /BLKR1/P I  ,  P I  ST,  P 1 1* 

COMMON  /BLKR4/RH00N , PH  I  ON . ZON , THE  I N , OMES , RNZ , VS , NF P 
COMMON  /BLKRA/RHOD . PH  I  0 , ZD , R I , STHE I N , CTHE I N 


.V.wv  'Vvr.-,' 
■ .  \ -vVv  ■ 


m 


FORMAT (2F7-1.3E10-3.5I5) 

FORMAT  (2F10.4,2E10.3,F10.4) 

FORMAT  (16 1  5) 

FORMAT  (El 0.3) 

FORMAT  (1  OF 7  - 1 ) 

FORMAT(8E10.3,/,4E10.3) 

FORMAT  (4E10.3) 

F0RMAT(19H  reference  temp.  =  , 1 1 X , F6 . 1 , 1 HK,/, 1 3H  temp,  span  =  , 

1 17X, F6. 1 , 1HK,/, 14H  capacitance  =,l6X,F10.5,  1 2HJ/cm/cm/cm/K , / , 

1  14H  conductance  =,l6x,F10.5.  6Hv/cm/K,/,  8Hlength  =,21X,F8.3, 

1 3H  cm,/,20H  #  of  nodal  points  =  , 1 1X , I  3. /» 1 3H  #  of  loops  =,l8X, 

1 1 3./.15H  #  of  unknows  =  , 1 6X , I  3 ,/. 1 2H  H  of  BC's  - . 1 9X . 13./. 

118h  #  of  temp.  BC's  =  . 1 3X, I  3 ,///) 

FORMAT (IX, 151 5) 

FORMAT (11  (1X.F5.1) ) 

F  ORMAT  (14(1X,F5*1)) 

FORMAT  (IX,  12F9.4) 

FORMAT  (IX, 12E9-2) 

WRITE (6,20) R£FT,REFDT,REFC,REFK,REFLC,NPT,NLP,NU,NBC,NBCT 

- <<< 

NORMALIZATION  «< 

«< 

- <<< 

REFQC=REFK*REFDT/REFLC 

REFA-REFK/REFC 

REFTK=REFLC*REFLC/REFA 

REFU=REFC*REFDT/REFTK 

WRITE  (6,11)  XC  (3)  , YC  (3)  ,C(3)  ,CK(3)  ,U0(3) 

DO  100  N-1.NPT 
C (N)  -C (N) /REFC 
CK  (N)-CK(N)  /REFK 
UO (N) - (UO (N) -REFT) /REFOT 
XC(N)-XC(N)  /REFLC 
YC(N)-YC(N)/REFLC 

WRITE  (6,11)  XC  (N)  ,  YC  (N)  ,C(N)  ,CK(N)  ,U0(N) 

- <<< 

«< 

MATRIX  FORMULATION  «< 

«< 

- <<< 

DO  300  1=1 , NLP 

WRITE  (6,21)  I  , LI  (I)  ,  II  (I)  ,J1  (I)  ,L2  (I)  ,  12  (I)  ,  J2  (I)  , 

1L3 (1)  .13(0  . J 3 ( 1 ) 

R1 2-DSQRT  ( (XC  (LI  ( I ) ) -XC (L2 ( I ) ) ) **2+  (YC  (LI  ( I ) ) -YC  (L2  ( I ) ) ) **2) 
R13-0SQRT((XC(L1  (I ) ) -XC  (L3  (D  ) )  **2+ (YC  (LI  (I) ) -YC  (L3  (I) ) )  **2) 
A12-(XC(L3(I))*(YC(L2(I))-YC(L1  (I) ) ) -YC  (L3  (D  )  *  (XC  (L2  (I) ) 

1-XC  (LI  ( I ) ) )  +XC  (L2  (I ) )  *YC  (LI  (I ) ) 

1  -XC  (LI  (l))*YC(L2(l)))*.500/R12 
A1 2-OABS  (A  1 2) 

A13*  (XC  (L2  (I) )  *  (YC  (L3  (I) )  "YC  (LI  (I) ) )  -  YC  (L2  (I) )  *  (XC  (L 3  ( I ) ) 

1-XC  (LI  (l)))+XC(L30))*YC(L1  (D) 

1  -XC  (LI  ( I ) )  *YC (L3 ( I ) ) ) * • 5D0/R1 3 
A1 3-DABS  (A13) 

V2*A1 2*R1 2* . 25D0 


V3-A13*R13*.25D0 


V1-A12*R12*.5D0 

CD1 2*  (CK  (LI  (l))+CK(L2(l)))*.5D0*A12/R12 

CD1 3=  (CK  (LI  ( I ) )  +CK  (L3(D))*.5D0*A13/R13 
IF  (  (11(1)  .GT.  0)  .AND.  (12(1)  .GT.  0)  .AND.  (13(1)  .GT.  0)  ) 
1  GOTO  200 

IF  (  (12(1)  .LE.  0)  .OR.  (13(1)  -LE.  0)  )  GOTO  120 

A (11  (I)  ,L1  ( I ) ) =A (LI  (I) ,L1  ( I ) ) -CD1 2-CD1 3 
B  (11  (I)  .  12  (l))=B  (LI  (I)  ,  12  (l))+CD12 
B  (LI  (I)  ,  I  3  ( I ) )  =B  (LI  (I)  ,  I3(I))+CD13 
D  (LI  (I) )  »D  (LI  (I)  )+C  (LI  (I) )  *V  1 
GOTO  200 

IF  (  (11(1)  .LE.  0)  .OR.  (I3(D  .LE .  0)  )  GOTO  130 

A  (L2  ( I )  ,  L2  (I ) )  =A  (L2  ( I )  .  L2  (I) )  -CD12 

B (L2 ( I )  ,11 ( I ) ) =B  (L2  ( I )  ,11  ( I ) ) +CD1 2 
D ( L 2 ( I ) )  =D (L2 ( I ) ) +C  (L2  ( I ) ) *V2 
GOTO  200 

IF  (  (11(1)  .LE.  0)  .OR.  (12(1)  .LE.  0)  )  GOTO  140 

A  (L3  ( I )  ,L3(I))-A(L3(I)  .  L3  ( I ) ) -CD1  3 
B  (L 3  ( I ) .11  ( I ) ) =B  (L3  (I ) .11 { I ) ) +CD1 3 
0(L3(I))-D(L3(I))+C(L3(D)*V3 
GOTO  200 

IF  (13(0  -LE.  0  )  GOTO  150 

A  (L 1  (I)  ,L1  ( I ) )  “A  (LI  (I)  ,L1  (I) ) -CD12-CD13 
A  (LI  (I)  ,L2(I))=A(L1  (I)  ,  L2  ( I ) )  +C01 2 
A  (L2  ( I )  ,L2(I))»A(L2(I)  , L2  ( I ) )  -CD1 2 
A  (L2  ( I )  ,L1  (l))-A(L2(l)  ,L1  (l))+CD12 
B  (LI  (I)  .  13(1))  =B  (LI  (I)  ,13(1))  +CD1 3 
D  (LI  (l))«D(L1  ( I ) )  +C  (LI  (l))*V1 
D(L2(I))-D(L2(I))+C(L2(I))*V2 
GOTO  200 

IF  (  12(1)  .LE.  0  )  GOTO  160 

A(L1  (I)  ,L1  (l))-A(L1  (I)  ,L1  (I) )  -CD12-CD1 3 
A  (LI  (I)  ,L30))«A(L1  (I)  , L 3  ( I ) )  +CD1  3 
A  (L 3  ( I )  .L3  (I) ) -A  (L3  (I)  .13  0))  -CD1 3 
A  (L3  ( I )  .LI  (I)) -A  (13(1)  ,L1  (I ) )  +CD1 3 
B  (LI  (I)  ,  I  2  ( I ) )  -B  (LI  (I)  ,12(1))  +CD1 2 
D  (LI  ( I ) )  *D  (LI  ( I ) )  +C  (LI  (l))*V1 
0(L3(I))-D(L3(I))+C(L3(I))*V3 

GOTO  200 

IF  (11(1)  .LE.  0  )  GOTO  170 

A  (L2  (I )  ,L2(I))-A(L2(I)  ,L2  (I) )  -CD12 
A (L 3 ( I )  * L3 ( I ) )  “A (L3 (I ) . L 3  ( I ) ) ~CD 1 3 
B  (L2  ( I )  ,11  (l))-B(L2(l)  .11  (DJ+CD12 
B  (L3  (I)  .11  (I) )  -B  (L3  (I)  .M  ( I ) )  +CD1  3 
D(L2(I))-D(L2(I))+C(L2(I))*V2 
D(L3(I))-D(L3(I))+C(L3(I))*V3 

GOTO  200 

A (LI  (I)  ,L1 (I ) ) *A  (LI  (I)  , L 1  (I) )  -CD12-CD13 
A  (LI  (I)  ,L2(I))-A(L1  (I)  ,  L2  ( I ) )  +CD1 2 
A  (LI  (I)  ,L3(I))-A(L1  (I)  ,  L 3  ( I ) )  +CD1  3 
A  (L2  ( I )  ,L2(I))-A(L2(I)  ,L2  (I) )  -C01  2 
A  (L2  (I)  ,  LI  (I) )  *A  (L2  (I)  ,  LI  (I ) )  +CD1 2 
A  (L3  ( I )  ,L3(I))-A(L3(I)  .  L3  ( I ) )  -CD1  3 
A  (L 3  ( I )  .LI  ( I ) ) "A  (L 3  ( I )  »L1  ( I ) ) +CD1 3 
0  (LI  (I) )  -D  (LI  ( I ) )  +C  (LI  (l))*V1 


D(L2(I))-0(L2(I))+C(L2(I))*V2 
D(L3(I))-D(L3(I))+C(L3(I))*V3 
200  IF  (  II  (I)  .GE.  0  )  GOTO  220 

11  (l)-IABS(ll  (I)) 

IF  (  (J1  ( I )  .NE.  L2  ( I ) )  .AND.  (J1(l)  .NE.  L2(I)+L3(D)  ) 

1  GOTO  210 

B  (LI  (I)  ,11  ( I ) )  =B  (LI  (I)  ,11  (l))+R12*.5D0 
210  IF  (  (J1(l)  .NE.  L3(U)  -AND.  (J 1  (!)  .NE.  L2(I)+L3(D)  ) 

1  GOTO  220 

B  (LI  (I)  ,11  (l))=B(L1  (I)  ,11  (l))+R13*.5DO 

220  IF  (  12  (I)  .GE.  0  )  GOTO  240 

12  (l)  =  IABS  (12(1)) 

IF  (  (J2(l)  .NE.  L  3  ( I ) )  .AND.  (J2  (I )  .NE.  L3  (D+L1  (I ) )  ) 

1  GOTO  230 

B(12(I),I2(I))=B(L2(I)  ,  12(1)) +R23*. 500 

230  IF  (  ( J  2  ( I )  .NE.  LI  ( I) )  .AND.  (J2(l)  -NE.  L3(D+L1(D)  ) 

1  GOTO  240 

B(L2(I)  ,  12(1))  =B(L2(I)  ,  12(1)) +R1  3*. 5D0 
240  IF  (  I3(D  -GE.  0  )  GOTO  260 

1 3  (I) -I  ABS  (13(1)) 

IF  (  ( J  3  (I)  .NE.  L1(D)  .AND.  (J3(D  .NE .  L1(I)+L2(D)  ) 

1  GOTO  250 

B(L3(D  .  1 3 ( D )  =B (L3 ( I )  ,I3(D)+R13*.5D0 
250  I  F  (  ( J  3  (I)  .NE.  L2  ( I ) )  -AND.  (J3(D  .NE.  L1(l)+L2(l))  ) 

1  GOTO  260 

B(L3(D  ,  13(D)  =B (L3 ( D  ,  13(1)) +R23*.5D0 
260  IF  (I  .LT.  NLP  )  GOTO  300 

WRITE  (6,22)  ((A  (K 1 , K2) ,K2=1 ,NU) ,K1=1 ,NU) 

WRITE  (6,23)  ( (B  (K 1 , K2)  ,K2=1,NBC)  ,K1«1,NU) 

300  CONTINUE 

DO  400  1*1, NU 
DO  320  J*1 , NU 

320  A  ( I  ,  J)  *A  ( I  ,  J)  /D  ( I) 

DO  3^0  J*1 , NBC 
340  B  ( I , J) *B  ( I , J) /D  ( I) 

400  CONTINUE 

WRITE  (6,22)  ( (A  (K 1 ,K2) .K2-1.NU) ,K1  =  1,NU) 

WRITE  (6,23)  ( (B  (K1 , K2) ,K2  =  1 ,NBC)  , K1=1 , NU) 

RETURN 

END 

c  ******  *  *  *  *  *  *  *  *  *  *  *  A  *  ********  A  A  *  *  A  *  *  *  *  *  it  *  *  *  *  * it  ;V  it  *  * it  A  it  ii  *  *  *  it  *  *  ******  *  *  *  * 
C***  *** 

C***  CONDUCTION  HEAT  FLUX  *** 

C***  *** 

c **************************************************** *************&**** 

SUBROUTINE  COND  (UB1 ,TKO,TK) 

IMPLICIT  REAL*8  (A-H ,0-2) 

COMMON  /BLKC2/A  (60,60) , B  (60, 60) 

COMMON  /BLKC3/D  (60)  ,  UBD  (60)  ,  UO  (60)  ,  UB  (60) 

COMMON  /BLKC4/REFT,REFDT,REFC,REFK,REFLC,NPT,NU,NBC,NBCT,NLP 
COMMON  /BLKC6/REFQC,REFA,REFTK,REFU 
D I  MENS  I  ON  UD  (60)  ,  UDO  (60)  ,  U  (60)  ,  UB1  (60) 

1  FORMAT  (IX, 1 2 F 9 •  4) 

2  FORMAT  (3F1 0.3) 


BOUNDARY  CONDITIONS 


C»> 

c»> - 

TK0=TKO/REFTK 
DO  1*20  1=1  ,  NBC 

IF  (l  .GT.  NBCT  )  GOTO  1*10 
UB (I)  =  (UB (I) -REFT) /REFDT 
GOTO  1*20 

1*10  UB  (l)-UB  (I)  /REFQC 

420  CONTINUE 

TK=TK/REFTK 
DTK=TK-TK0 
TK0=TK 

DO  530  1=1 ,NBC 

IF  (I  .GT.  NBCT  )  GOTO  510 
UB1  ( I )  =  (UB1  (I) -REFT) /REFDT 
GOTO  520 

510  UB1  (I)  =UB1  (I)  /REFQC 

520  UBD(I)*UB1  (l)-UB(l) 

530  UB(I)-UBKI) 

WRITE  (6,1)  (UB1  (I) , 1=1 , NBC) 

WRITE (6 , 2) TKO.TK , DTK 

C»> - - - - - 

c»> 

C»>  SERIES  EXPANSION 

C»> 

c»> - 

DO  620  1=1, NU 
UDO  (l)-O.DO 
DO  610  J*1 , NU 

610  UDO  (I)  -UDO  (l)-t-A  (I  ,J)*UO(J) 

DO  620  J-1.NBC 

620  UOO  (I)  =UDO  (l)+B  (I  ,J)  *UB  (J) 

DO  630  1=1, NU 

UDO  ( I ) =UDO  ( I ) *DTK 

630.  U  (I)  =U0  (l)+UDO  (I) 

WRITE  (6,1)  (U  ( I )  ,1*1, NU) 

DO  650  1=1 ,NU 
UD  ( I ) =0 . DO 
DO  640  J*1 , NU 

640  UD(l)-UD(l)+A(l,J)*UDO(J) 

DO  650  J=1 , NBC 

650  UD  (l)»UD  (l)+B  (I  ,  J)  *UBD  (J) 

DO  700  1=2,50 
TEST-O.DO 
DO  660  J=1 ,NU 

UD  (J)»UO  (J)  ftDTK/l 
TEST-TEST+UO  (J)*UD(J) 

UDO (J) =UD  (J) 

660  U  (J)  =U  (J)  +UD  (J) 

WRITE  (6,1)  (U(J)  , J=1 , NU) 
TEST-OSQRT  (TEST/NU) 

IF  (  TEST  .LE.  .10-10  )  GOTO  800 
DO  670  K=1 , NU 


0.0000  0.0000  .350E+001  .392E+000  500.0000 

0.5000  0.0000  -350E+001  .392E+000  500.0000 

1.0000  0.0000  . 350E+001  .392E+000  500.0000 

1.5000  0.0000  .350E+001  .392E+000  500,0000 

2.0000  0.0000  .350E+001  .392E+000  500.0000 

2.5000  0.0000  -350E+001  .392E+000  500.0000 

3.0000  0.0000  .350E+001  . 392E+000  500.0000 

3.5000  0.0000  -350E+001  .392E+000  500.0000 

it. 0000  0.0000  .350E+001  .392E+OOO  500.0000 

it. 5000  0.0000  -350E+001  . 392E+000  500.0000 

5.0000  0.0000  . 1 OOE+OOI  . 392E+000  500.0000 

1-15  7  5  '16  1  2-26  1 


at 


UD  (K) =0 .DO 
DO  670  L=1 ,NU 
UD  (K)  =UD  (K)  +A  (K ,  L)  *UD0  (L) 
CONTINUE 
DO  810  1=1, NU 
UO  (I)  =U  (1) 


RETURN 

END 

//GO. SYS  IN  DO  * 

1000.0  4. 921 


22.326 


500.0  500.0  . 1 00E+001  . 1 00E+001  .500E+001 


0.0000 
0.0000 
0.0000 
0.0000 
0.5000 
0.5000 
0.5000 
0.5000 
1 .0000 
1 .0000 
1.0000 
1 .0000 
1 .5000 
1 .5000 
1 .5000 

1 .5000 
2.0000 
2.0000 
2.0000 
2.0000 

2.5000 
2.5000 
2.5000 

2.5000 
3.0000 
3.0000 
3.0000 
3.0000 

3.5000 
3.5000 
3.5000 

3.5000 

4.0000 

4.0000 

4.0000 

4.0000 

4.5000 
4.5000 
4.5000 
4.5000 
5.0000 
5.0000 
5.0000 
5.0000 


2.0000  .350E+001  .392E+000  500. 

1.5000  .350E+001  .392E+000  500. 

1.0000  .350E+001  .392E+000  500. 

0.5000  .350E+001  .392E+000  500. 

2.0000  .350E+001  .392E+000  500. 

1.5000  .350E+001  . 392E+000  500. 

1.0000  .350E+001  .392E+000  500. 

0.5000  .350E+001  .392E+000  500. 

2.0000  .350E+001  .392E+000  500. 

1.5000  .350E+001  . 392E+000  500. 

1.0000  .350E+001  . 392E+000  500. 

0.5000  .350E+001  .392E+000  500. 


2.0000  .350E+001  . 392E+000 
1.5000  .350E+001  .392E+000 
1.0000  .350E+001  .392E+000 
0.5000  .350E+001  .392E+OOO 
2.0000  .350E+001  .392E-t-000 

1.5000  .350E+001  .392E+000 


.392E+000 

■392E+000 


■392E+000  •  500. 


1.0000  .350E+001  .392E+000  500. 
0.5000  .350E+001  .392E+000  500. 


2.0000  .350E+001 
1.5000  .350E+001 


•392E+000  500. 
.392E+000  500. 


1.0000  .350E+001  .392E+000  500. 
0.5000  .35OE+OOI  .392E+000  500. 


2.0000  .35OE+OOI  .392E+000 
1.5000  .350E+001  . 392E+000 
1.0000  .350E+001  . 392E+000 
0.5000  .350E+001  . 392E+000 
2.0000  .350E+001  .392E+000 
1.5000  .350E+001  . 392E+000 
1.0000  .350E+001  .392E+000 
0.5000  .350E+001  . 392E+000 
2.0000  .350E+001  .392E+000 
1.5000  .350E+001  .392E+000 

1.0000  .350E+001  .392E+000 

0.5000  .350E+001  . 392E+000 
2.0000  .350E+001  .392E+000 
1.5000  .350E+001  .392E+000 
1.0000  .350E+001  . 392E+000 
0.5000  .350E+001  . 392E+000 
2.0000  .350E+001  .392E+000 
1.5000  -350E+001  . 392E+000 

1.0000  .350E+001  .392E+000 

0.5000  .350E+001  . 392E+000 


1 .0000 

0.5000 


392E+000  500. 

392E+000  500. 

392E+000  500. 

392E+000  500. 

392E+000  500. 
392E+000  500. 

.392E+000  500. 
•392E+000  500. 

392E+000  500. 
392E+000  500. 
392E+000  500. 

392E+000  500. 

392E+000  500. 
392E+000  500. 

392E+000  500. 

392E+000  500. 

392E+000  500. 

392E+000  500. 


16.0 
55  80 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000' 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 
0000 


34.0  10004112.0 
44  28  11 


v'v  i-':- 


•  A  • 


•500E+03  .500E+03  .500E+03  . OOOE+OO  . OOOE+OO  .OOOE+OO  .202E+00  .174E+01 
.2&7E+01  . 344E+01  .409E+01  .486E+01  .522E+01  -575E+01  .623E+01  .671E+01 
.360E+01  .OOOE+OO  .OOOE+OO  .OOOE+OO 

1000.0000  0003 .4417  0023.0000  0004.9213  0005.7150  0022.3266 
0000.8000  5.080E-03  0015.0000  0000.0000  0016  0034  0000  0004 

0002  0002  0002.5400  0000. OOCO  0000.0000 

0500.0000  0500.0000  0500.0000  0500.0000  0000.9500  1 . 0OOE+O4 

0170.0000 

0052.4007  0.800E-05  0007.6500  0075-6400  0038.0000 

0001.0000  0007.6500  0075.6400  0038.0000 

0000.0100  0007.6500  0075.6400  0038.0000 

0000.0100  0007.6500  0075.6400  0038.0000 

.1000 

.500E+03  .500E+03  .500E+03  .500E+03  .500E+03  .500E+03  .500E+03  .500E+03 

.500E+03  .500E+03  .5OOE+03  .OOOE+OO  .OOOE+OO  .OOOE+OO  .202E+00  . 1 74E+01 

. 267E+01  . 344E+01  .409E+01  .486E+01  . 522  E+01  .575E+OI  .623E+01  .671E+01 

•360E+01  .OOOE+OO  .OOOE+OO  .OOOE+OO 

1000.0000  0003.4417  0023.0000  0004.9213  0005.7150  0022.3266 
0000.8000  5.080E-03  0015.0000  0000.0000  0016  0034  0000  0004 
0002  0002  0002.5400  0000.0000  0000.0000 

0500.0000  0500.0000  0500.0000  0500.0000  0000.9500  1 . 000E+04 
0180.0000 

0078.7007  0.800E-05  0005.3000  0051.2800  0026.0000 

0001.0000  0005.3000  0051.2800  0026.0000 

0000.0100  0005.3000  0051.2800  0026.0000 

0000.0100  0005.3000  0051.2800  0026.0000 

.1000 

.500E+03  .500E+03  .500E+03  .500E+03  .500E+03  .500E+03  .500E+03  .500E+03 

.500E+03  .500E+03  .500E+03  .OOOE+OO  .OOOE+OO  .OOOE+OO  .202E+00  .174E+01 

.267E+01  . 344E+01  .409E+01  .486E+01  .522E+01  .575E+01  .623E+OI  .671 E+01 

•360E+01  .OOOE+OO  .OOOE+OO  .OOOE+OO 

1000.0000  0003.4417  0023.0000  0004.9213  0005.7150  0022.3266 
0000.8000  5.O8OE-03  0015.0000  0000.0000  0016  0034  0000  0004 
0002  0002  0002.5400  0000.0000  0000.0000 
0500.0000  0500.0000  0500.0000  0500.0000  0000.9500  1 . OOOE+04 
0190.0000 

0085.6907  0.800E-05  0002.9500  0026.9200  0014.0000 

0001.0000  0002.9500  0026.9200  0014.0000 

0000.0100  0002.9500  0026.9200  0014.0000 

0000.0100  0002.9500  0026.9200  0014.0000 

.1000 

.500E+03  .500E+03  .500E+03  .500E+03  .500E+03  .500E+03  .500E+03  .500E+03 

.500E+03  .500E+03  .5OOE+O3  .OOOE+OO  .OOOE+OO  .OOOE+OO  .202E+00  . 1 74E+01 

. 267E+01  . 344E+01  .409E+01  .486E+01  . 522E+01  -575E+01  .623E+OI  .671E+01 

•360E+01  .OOOE+OO  .OOOE+OO  .OOOE+OO 

1000.0000  0003.4417  0023.0000  0004.9213  0005.7150  0022.3266 
0000.8000  5.080E-03  0015.0000  0000.0000  0016  0034  0000  0004 
0002  0002  0002.5400  0000.0000  0000.0000 
0500.0000  0500.0000  0500.0000  0500.0000  0000.9500  1. 000E+04 
0200.0000 

0073-8000  0.800E-05  0000.6000  0002.5600  0002.0000 

0001.0000  0000.6000  0002.5600  0002.0000 

0000.0100  0000.6000  0002.5600  0002.0000 

0000.0100  0000.6000  0002.5600  0002.0000 

.1000 


52 


34 

0 

0 

30 

0 

0 

35 

31 

0 

0 

30 

0 

0 

35 

31 

0 

0 

35 

0 

0 

32 

36 

0 

0 

32 

0 

0 

35 

37 

-24 

33 

33 

-23 

37 

38 

34 

0 

0 

33 

-23 

0 

38 

34 

0 

0 

35 

0 

0 

38 

39 

0 

0 

35 

0 

0 

38 

39 

0 

0 

35 

0 

0 

40 

36 

0 

0 

35 

0 

0 

40 

41 

-25 

79 

42 

-12 

41 

38 

0 

0 

42 

0 

0 

42 

-12 

43 

43 

-13 

42 

38 

39 

0 

0 

43 

-13 

0 

38 

43 

-13 

44 

44 

-14 

43 

39 

40 

0 

0 

39 

0 

0 

44 

4 

-28 

45 

45 

1 

0 

8 

46 

2 

0 

45 

1 

0 

8 

46 

2 

0 

47 

3 

0 

8 

12 

0 

0 

8 

0 

0 

47 

47 

3 

0 

12 

0 

0 

48 

16 

0 

0 

12 

0 

0 

48 

20 

0 

0 

16 

0 

0 

49 

48 

4 

0 

49 

5 

0 

16 

24 

0 

0 

20 

0 

0 

50 

49 

5 

0 

20 

0 

0 

50 

28 

0 

0 

24 

0 

0 

51 

50 

6 

0 

24 

0 

0 

51 

28 

0 

0 

32 

0 

0 

51 

52 

8 

0 

51 

7 

0 

32 

36 

0 

0 

53 

9 

0 

32 

52 

8 

0 

53 

9 

0 

32 

40 

0 

0 

54 

10 

0 

36 

53 

9 

0 

54 

10 

0 

36 

44 

-14 

55 

40 

0 

0 

55 

54 

10 

0 

55 

1 1 

0 

40 

151 .000  1 .000  50  92484 


101.325  300.000  50.0  151.000  i.ooo 

1 

1715.61  1805.19  1901.64  1991.21  2094.56  2273' 
3026.46  3217.63  3410.55  3610.36  3844.62  3851. 
5408.65  5684.25  5959.85  6235.45  6511.05  6779. 
7928.10  8095.75  8233.55  8302.45  8350.68  8371 ' 
8398.91  8392.02  8378.24  8343.79  8260.00  8164, 
1000.0000  0003.4417  0023.0000  0004.9213  0005, 
0000.8000  5.080E-03  0015.0000  0000.0000  0016 
0002  0002  0002.5400  0000.0000  0000.0000 

0500.0000  0500.0000  0500.0000  0500.0000  0000, 
0160.0000 

0029.1007  0.800E-05  0010.0000  0100.0000  0050, 
0001.0000  0010.0000  0100.0000  0050, 
0000.0100  0010.0000  0100.0000  0050, 
0000.0100  0010.0000  0100.0000  0050, 

.1000 

.500E+03  .500E+03  .500E+03  .500E+03  .500E+03 


70  2389.94  2501.07  2631.98  2859.35 
00  4092.66  4347.59  4878.12  5139.94 
76  7027.80  7068.95  7510.10  7716.00 
35  8378.24  8385.13  8392.02  8398.91 
65  8026.85  7854.60  7579.00  7234.00 
7150  0022.3266 
0034  0000  0004 

9500  1  . 000E+04 


.500E+03  .500E+03  .500E+03 


.500E+03  .500E+03 
.500E+03  .5O0E+O3 

.267E+01  .344E+01 
.36OE+OI  . OOOE+OO 

/* 


•500E+03  .500E+03 
.500E+03  .OOOE+OO 
.  409E+01  .486E+01 
.OOOE+OO  .OOOE+OO 


.500E+03  .500E+03 
.OOOE+OO  .OOOE+OO 
.522  E+ 0 1  . 575E+01 


.500E+03  .500E+03 

.  202E+00  .174E+01 
.623E+01  .671E+01 


53 


